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Abstract 


An approach for modeling and simulating landing gear systems is presented. 
Specifically, a nonlinear model of an A-6 Intruder Main Gear is developed, simulated, 
and validated against static and dynamic test data. This model includes nonlinear effects 
such as a polytropic gas law, velocity squared damping, a geometry governed model for 
the discharge coefficients, stick-slip friction effects and a nonlinear tire spring and 
damping model. An Adams-Moulton predictor corrector was used to integrate the 
equations of motion until a discontinuity caused by a stick-slip friction model was 
reached, at which point, a Runga-Kutta routine integrated past the discontinuity and 
returned the problem solution back to the predictor corrector. Run times of this software 
are around 2 minutes per 1 second of simulation under dynamic circumstances. To 
validate the model, engineers at the Aircraft Landing Dynamics facilities at NASA 
Langley Research Center installed one A-6 main gear on a drop carriage and used a 
hydraulic shaker table to provide simulated runway inputs to the gear. Model parameters 
were tuned to match one dynamic case. Other cases were then run with the updated 
parameters and the results were in excellent agreement with the test data. 
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Chapter 1: Introduction 


1.1 Background 

In recent years, NASA and many aerospace companies have increased their 
research focus on the High Speed Civil Transport (HSCT). The concept is to fly a 
supersonic (mach 2.4) airplane to various places on the globe at an economical price for 
both carrier and passenger use. Its overall appearance will be similar to that of the 
Concorde, a current, expensive supersonic carrier. 

To make the HSCT more cost effective, much effort is being expended in the 
design stage. One major problem encountered in its development is the trade-off between 
structural rigidity and total weight. To this point, the structure of the fuselage and wings 
has been designed for aerodynamic performance. In the early stages of design, however, 
landing gear location and dynamic performance are rarely considered. Since the fuselage 
on the HSCT is very long and slender, it is very sensitive to external, low frequency 
disturbances, or vibrations. Therefore, a goal of the landing gear design is to reduce 
disturbance transmission from the ground to the fuselage. 

Computer simulations are being developed to study this disturbance transmission 
problem. The task is to take information concerning gear dynamics, fuselage dynamics, 
runway profile, and taxi speed, to develop a simulator for predicting aircraft ground 
response. Simulation of landing gear dynamics has been the subject of much research for 
many years. The military has long been interested in simulating gear response to 
repaired, bomb-damaged runways. A great deal of effort has been applied to the problem 
of determining how well to repair a runway to prevent landing gear failures. This effort 
did not focus on changing the gear (i.e. to control the force transmission), but rather on 
changing the runway repair specifications. Active control concepts may render landing 
gears less sensitive to rough runways, decreasing the time needed to repair damaged 
runways, and thus allowing quicker response of military missions. 
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1.2 Objective 

This document will present an approach for modeling and simulating landing gear 
systems. Specifically, a nonlinear model of an A-6 Intruder Main Gear is developed, 
simulated, and validated against test data. This model includes nonlinear effects such as a 
polytropic gas model, velocity squared damping, a geometry governed model for the 
discharge coefficients, stick-slip friction effects and a nonlinear tire spring and damping 
model. To validate the model, engineers at the Aircraft Landing Dynamics Facility at 
NASA Langley Research Center installed one A-6 main gear on a drop carriage and used 
a hydraulic shaker table to provide simulated runway inputs to the gear. Model validation 
used both quasi-static and dynamic tests. In summary, then, this research presents a 
comprehensive mathematical formulation of landing gear systems, verifies the modeling 
techniques with tests, and discusses approaches for further model correlation using the 
test results. 

13 Literature Survey 

Concurrent and past work of this kind have been generally to predict taxi loads of 
military aircraft over repaired, bomb-damaged runways. This work has led to extensive 
modeling of military aircraft with the goal of determining minimum repair procedures to 
runways to prevent gear failure on rollout. One major accomplishment of this research is 
the HAVE BOUNCE 1 simulation program. Using this program, the USAF has 
simulated the dynamic response of many military aircraft over bomb damaged runways. 
These simulations are validated with test data and are used to identify component 
weaknesses and operational limits. Validation of these simulations is usually achieved 
through the use of the Aircraft Ground Induced Loads Excitation (AGILE) 2 test facility at 
Wright-Patterson Air Force Base. Simulations usually combine nonlinear coupled 
differential equations of the landing gear with linear techniques describing the fuselage 
structure to generate the total aircraft response. Each of these simulations is usually very 
good, but each is also very tailored for the plane in question. In addition, the militaiy is 
concerned mainly with the problem of traversing repaired sections sequentially on a 
runway. They have, therefore, limited the inputs to their models mainly to the various 
classes of repairs. 
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Another approach discussed in the literature uses data from tests to determine 
parameters of a state space model for the landing gear and fuselage 3 . These parameters 
are included in a quasi-linear formulation through look-up tables. Depending on the exact 
form of the model and the measurement capability, the results of this type of formulation 
are fair to good in comparison to actual test data. This approach is good when test data is 
easily accessible for parameter determination, and simulation computation time is limited. 
Freymann 4 revisited an experiment by Ross and Edson 5 in which an actively controlled 
servovalve is connected to a landing gear to augment damping in the system. Ross and 
Edson 5 described an electronic controller and an actively controlled landing gear which 
was found to significantly reduce forces sustained by an aircraft during takeoff, landing 
impact and rollout. The results were obtained analytically through the use of a linearized 
model of the equations of motion and confirmed experimentally. The servo-controller 
was designed to maintain a certain command force level by porting flow into and out of 
the landing gear. This approach to active control, however, is very expensive in terms of 
weight and complexity. There, also, was no recourse developed in case of servo failure. 
Their model included a linear tire spring, no tire damping, no metering pin and no 
rebound chamber. However, they did include friction and velocity squared fluid 
damping. Freymann extended this research by implementing this active control scheme 
on a fighter aircraft and testing this aircraft's response to various frequency inputs. The 
control system performed well in attenuating aircraft motions. 

Much effort has been exerted by Stirling Dynamics 6,7 in the field of simulation 
and control. Their main simulation program includes many highly detailed models of the 
various components of several types of landing gear. The program allows part selection 
for individual gear types and outputs landing gear dynamic responses in terms of 
positions, velocities, accelerations and angular equivalents for any part of the gear or 
fuselage which is selected as a node. Other outputs were subsystem force interactions. 
Extensive validation has proven this to be a very accurate simulation tool. The complexity 
of the model generally captured most dynamic effects seen in test data. This same 
program was used to test an active orifice concept applied to the nose gear of a typical 
transport plane. These results showed reduction of peak and root mean square values of 
normal acceleration, especially in the nose and tail section of the plane. An active orifice 
mechanism was not detailed, only an assumed behavior. This is a very good program for 
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landing gear design and testing of existing gear configurations. However, a simpler 
model would allow the physics of the strut, only, to be scrutinized, leading perhaps to a 
clearer understanding of landing gear behavior. 

Research into the behavior of a supersonic carrier during ground operations was 
performed by C. G. Mitchell 8 . His theoretical analysis and test experience with the 
Concorde has shown that the supersonic transport is more sensitive to uneven runways 
than the subsonic transport. Results reported in [8] show that much care must be taken to 
minimize undercarriage stiffness and friction if problems of cockpit vibrations and 
airframe and undercarriage fatigue were to be avoided. 

Ramamoorthy 9 performed a parameter study on his model of an articulated 
landing gear and found that changes in the discharge coefficient could alter the results 
dramatically. No quantitative conclusions were made and no validation of the model was 
performed. However, Wahi 10 found that the coefficient can alter forces transmitted to the 
fuselage by as much as 25% and that proper estimation of this parameter needs to be 
based on both the Reynolds number and orifice geometry. A semi-empirical model 
developed by Bell, Schlichting, Knudsen et. al. was used for this parameter. With this 
model for the discharge coefficient, Wahi 10 developed a landing gear model with two 
degrees of freedom to investigate the optimization of the metering pin shape. The results 
of this model compared well to flight and drop tests. Optimization of the metering pin, 
for this particular case, showed some improvement in the force reduction. Once a 
metering pin shape has been defined, it cannot change. An active orifice concept would 
allow the damping characteristics of the gear to be continuously changed in reaction to 
any input to reduce vibration transmission. 

An optimization of many strut characteristics was performed by Li, Gou-zhu, and 
Qing-zhi”. This paper described an optimization approach for landing gear design using 
as design variables the initial pressure, initial air volume, and an artificial oil damping 
coefficient (which in reality, is a function of the hydraulic and pneumatic areas as well as 
the discharge coefficient and oil density). The objective function was the mean square 
value of the fatigue power spectral density. Constraints were in the form of landing 
impact energy, static compression ratio, maximum compression ratio and limits on the 
damping coefficient. The results show a significant reduction in the accumulated fatigue 
on the simulated Boeing 707. What the results fail to show, however, is whether the 
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upper and lower limits on the damping coefficients are physically achievable. The limits 
did not include considerations of geometry and realistic discharge coefficient values. 

Doyle 12 provides an excellent literature survey on aircraft ground dynamic 
simulation techniques. His report contains a brief summary of the computer programs 
written to predict the dynamic displacements and forces resulting from nonflight aircraft 
operations. The capabilities of each program and their limitations and numerical 
techniques are cited. 


1.4 Research Significance 

The significance of the material treated in this research is that it brings together in 
one place a comprehensive development of the theory of telescoping gear. This 
document contains the development of the equations of motion and details the more 
standard practices of expressing them in terms of physically measurable quantities. The 
model has only two degrees of freedom, both in the vertical direction. In the investigation 
of loads transmitted into the fuselage, though, this is the most important direction. The 
model is fully nonlinear and includes such effects as a polytropic gas model, a velocity 
squared damping term, which includes a discharge coefficient that is a function of orifice 
geometry, extension damping, stick-slip friction in the gear, and nonlinear tire model. All 
parameters such as polytropic gas constant, orifice geometry, frictional quantities, etc. 
appear explicitly in the equations, and can be used in a sensitivity analysis. Also, 
optimization of gear geometry and initial charge pressures and volumes is easily 
accomplished using this model. In the end, control concepts can be linked to this model 
for investigation of force transmission reduction. 

This research also treats the subject of numerical integration of the equations of 
motion. The stiff, nonlinear, and discontinuous behavior of these equations make this a 
difficult problem to solve numerically, and many considerations were made to make it 
easier. Also, this document details a series of tests and procedures by which to validate 
the model. This validation is both static and dynamic. A frequency response method 
was used to update the parameters in the model and other types of cases to validate the 
simulation. 
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1.5 Document Outline 

This document is divided into six chapters. After the introduction in Chapter 1, 
Chapter 2, discusses the theoretical and mathematical development of the equations of 
motion of a landing gear. Chapter 3 details the method in which these equations were 
numerically implemented and some of the problems that were encountered. Chapter 4 
discusses the equipment used in the experimental validation effort. The next chapter 
describes the test procedures that were implemented to validate this simulation and some 
statements about error control. Finally, Chapter 6 will discuss some future research and 
experimental plans and presents some concluding remarks. Also included, in Appendix 
A, is the FORTRAN program used to obtain the simulation results shown in this 
document. 
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Chapter 2: Problem Formulation 


2.1 Initial Landing Gear Investigation 

This chapter is intended to familiarize the reader with landing gear terminology 
and to demonstrate a mathematical development of the equations of motion for a 
telescoping landing gear. Figure 2-1 is intended to acquaint the reader with basic landing 
gear components. It shows the simplified components of a telescoping, main landing 
gear (as opposed to a nose gear). 



1) Upper Mass (Fuselage) 

2) Nitrogen Gas (Pneumatic Spring) 

3) Outer Cylinder 

4) Hydraulic Fluid 

5) Orifice Plate 

6) Metering Pin 

7) Snubber Orifice 

8) Snubber (Rebound) Chamber 

9) Lower Piston 

10) Tire 


Figure 2-1: Schematic of typical telescoping main landing gear. 
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Point 1 on the figure is a rigid body representation of the aircraft fuselage that the 
gear carries and is the interface between the plane and the gear. Point 2 is a chamber 
containing compressed nitrogen which serves as a spring that carries the weight of the 
plane in ground operations. Point 3 refers to the main, upper cylinder which houses the 
compressed gas, hydraulic fluid, and within which the piston slides. The hydraulic fluid 
is represented by the shaded area located by point 4. Point 5 is the orifice plate. It is 
essentially a circular plate with a hole in the center through which the hydraulic fluid 
flows when the strut is stroking. It, along with the metering pin, point 6, controls the 
damping characteristics of the gear. The metering pin is rigidly fixed to the piston head. 
As the strut strokes, the changing size of the metering pin passes through the constant 
hole in the orifice plate, causing a variable effective orifice diameter, i.e. variable fluid 
damping. Nose gear on most planes have no metering pin. Point 7 locates one of many 
rebound or snubber orifices (usually around 12, depending on the gear). These holes lead 
into a small volume on the backside of the piston head (point 8) called the rebound or 
snubber chamber. The purpose of the snubber is to provide damping when the strut 
extends. The snubber orifices are variable in that they are dependent upon a slip ring that 
either allows a large orifice in the compression stage or a smaller orifice in the extension 
stage. Point 9 is the piston. It houses the metering pin and is also the rigid connection of 
the wheel axle. Finally, point 10 is the tire. This element of the gear adds both spring and 
damping characteristics to the overall performance of the gear, and is selected carefully 
for various applications. 

A study by Ross and Edson 5 of Hydraulic Research provided the initial 
information for this research. They developed nonlinear equations for a simplified 
telescoping landing gear. They then linearized the equations about the ground equilibrium 
point and studied the effect of an active damping control scheme. The model they 
developed did not include a metering pin or a rebound chamber. However, their model 
included a servovalve to port fluid from one chamber (upper or lower) to the other to 
control the gear damping response, precluding the usefulness of the snubber or metering 
pin. This method of control requires high hydraulic pressures and large pumps and 
plumbing to accomplish the task, making it difficult to implement. The report, however, 
did conclude that active control gear can reduce the ground loads transmitted to the 
fuselage. Their simulation was validated using experimental equipment and facilities of 
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Hydraulic Research. The linearized model did not allow explicit investigation of orifice 
diameters and other parameters critical in understanding the landing gear's dynamic 
behavior. But, these studies were considered investigative in nature and led to a more 
complete understanding of some of the complex dynamics of landing gear. 

2.2 Nonlinear Model Development 

To extend the work by Ross and Edson 5 , this research discusses an independent 
development of a mathematical model of a main landing gear with all the relevant 
physical parameters included. The nonlinear equations of motion are developed for a 
telescoping main gear. The analytical model used is a representation of an A-6 Intruder 
main gear. This gear was chosen because facilities exist to test and characterize an A-6 
gear for simulation validation. Specific details of the gear were taken from the technical 
drawings supplied by the Grumman Company. 

An initial model was developed that only included the air-spring above the fluid, 
fluid dynamics through a fixed orifice, and a linear tire spring term. This simple model 
allowed some trend comparison between the results of this model and the early results of 
the linearized gear of Edson and Ross 5 . A metering pin was then added to change the 
main orifice effective diameter as a function of stroke. Another variation from Edson and 
Ross was the addition of a snubber, or rebound chamber. This feature provides damping 
while the gear is extending. Since this new model is to be validated with test data, some 
attempt to quantify frictional effects was also made. The model includes constant seal 
friction as well as a variable friction that is a function of stroke. In a further effort to be 
realistic, a nonlinear tire model was added. This tire model has a spring rate that is a 
function of tire deflection and damping proportional to compression rate. In the equations 
developed below, the spring and damping coefficient are used as if they were constant. 
The nonlinear characteristics of each of these terms is included in the equations of motion 
that are actually integrated. 

Figure 2-2 is a schematic of the gear used in the development of the equations of 
motion. This schematic is representative of a general telescoping-type main landing gear. 
It includes the aerodynamic lift on the plane, Lift, the upper mass (of the plane's fuselage) 
and the mass of the main cylinder lumped together as a rigid mass, and the mass of 
the piston and the mass of the tire, also lumped together as M L . The inertial coordinate of 
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the upper mass is X wg . The zero value for X wg is when the gear is fully extended with the 
tire just touching the ground. From this same gear configuration, X a , the coordinate of 
the lower mass, is taken as zero at the axle of the tire. Therefore, when the gear is in 
some compressed state, X, measures the deflection of the tire when the ground input, 

U(t), is zero. In the compressed nitrogen chamber (upper cylinder) with cross sectional 
area of Ay, the pressure is P u . Likewise, in the lower chamber with cross sectional area of 
A l , there is a pressure of P L . In the snubber chamber, with annulus area of A R , the 
pressure is defined to be P s . The orifice plate has a hole of diameter through which 

the metering pin, with variable diameter D pin moves. Fluid reaches the snubber chamber 
through the orifices d s c and d s E , where the superscripts represent either the compression 
mode or extension mode respectively. The diameter of the piston, D pi , is used to calculate 
A r . Simply subtract the area of the piston shaft from that of the lower cylinder to get A„. 
The tire is also shown in Figure 2-2 with a distinction of pointing out that the tire spring 
and damping coefficients, K,, and C, are nonlinear and contribute to the calculation of the 
tire force F,. 


10 



Upper Mass/ Cylinder M 



Slip Ring Snubber 
d s c Large 
(Compression) 
d s E Small 
(Extension) 


+X 


"g 


Piston/ Lower 
Mass 


■ ; J; - ' 

■ ■ ■■■• 

s 


K„ C t =>F, 
Nonlinear 


+U(t) t 



Figure 2-2: Schematic of a telescoping main landing gear. 
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Figure 2-3: Schematic of upper mass and main cylinder 

Figure 2-3 shows the forces acting on the upper mass. Balancing the forces on the upper 
mass gives the following equation: 

=Kg-L- PAo - P l{\ - A a) +P s\ + f ( 2 . 1 ) 

The term on the left hand side of Eq. (2.1) is the inertial motion term, g is the gravitational 
acceleration, f is the friction present in the gear, and all other terms are as described 
previously. This equation assumes that the fluid pressure in the upper cylinder is identical 
to the pneumatic pressure. In this development, the variable A c , the main orifice area, 
reflects the fact that the metering pin is included, i.e. it is a variable cross-sectional area 
depending on stroke. 
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Figure 2-4: Schematic of lower mass. 

Figure 2-4 shows the forces acting on the piston. Summing the forces on the lower mass 
(piston) the force balance equation is: 

M ( X a = M,g + P L (A L - A s )~ P s (A r -A s )~ F'± f (2.2) 

where the left hand side of Eq. (2.2) is the inertial motion of the lower mass and A s is the 
area of the snubber orifice. F, is the force that is transmitted through the tire from the 
ground and has the form: 

F, = K,(x a +u)+c,(x„ + u) 

where the tire force is a function of a nonlinear tire stiffness and a damping force that is 
composed of a damping coefficient that is proportional to the tire stiffness and the time rate 
of change of the tire deflection. 

2.3 Relation of Pressures to Stroke Position and Stroke Rate 

The pressure terms in Eqs. (2.1) and (2.2) are as yet unknown and need to be 
related to the positional variables X wg and X a or their derivatives. The pressure of the 
compressed nitrogen in the upper cylinder can be described by the polytropic gas law for a 
closed system as: 
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(2.3) 



where X s is the stroke available, given by: 

X s =X wg - X a (2.4) 

with X S1 as some initial length, P^, the charge pressure at X^, and y, the polytropic gas 
constant. X smax is the maximum value to which the gear can be extended. This form of 
representation of the pressure change is assumed to happen as a quasi-equilibrium 
process 13 . The significance of the polytropic gas constant is that it describes the type of 
process that occurs. For example, if y = 0, the process would be isobaric, or constant 
pressure. However, for an ideal gas, y = 1 , corresponding to an isothermal (constant- 
temperature) process. If y is equal to the ratio of the specific heats of a gas, the process is 
isentropic, or constant entropy. These cases are idealizations of particular processes. In real 
situations though, the polytropic gas constant is not constant at all and is usually calculated 
from pressure-stroke data. An average value is usually sufficient in application. 

Equation (2.3) was defined in such a manner that P u will become very large when 
X s is near X smju , i.e. the gear is nearly completely collapsed. This is a suitable 
representation of the process, with only the polytropic gas constant y as an unknown. 

The pressures (P L and P s ) of the fluid in the lower cylinder and in the snubber are 
related to the flow rates of the fluid into and out of those regions. The volumetric flow 
rates through the orifice plate hole, Q c , and the snubber orifices, Q s , can be determined by 
combining the continuity equation and Bernoulli's equation for fluids 14 . Flow is always 
from the higher pressure to the lower pressure. Bernoulli's equation for an incompressible 
fluid states that along a streamline 14 , 

P/d + ( l/2g)V 2 + Z = Constant (2.5) 

where P is the pressure at some point, g is the gravitational acceleration, V is the velocity of 
the flow, d is the specific weight of the fluid which is equal to the fluid density (p) 
multiplied by the gravitational acceleration (g), and Z is the height difference from some 
zero reference. This equation assumes that the viscous effects within the fluid are 
negligible, the flow to be steady and incompressible, and that the equation is applicable 
along a streamline. 

Equating Bernoulli's equation (Eq. (2.5)) at two points in the flow along the same 
streamline yields: 


14 



P,A) + (l/2g)V, 2 + Z, = P/o + (l/2g)V 2 2 + Z, (2.6) 

In the case of a landing gear, the potential distance between Z, and 2^ can be neglected as 
the distances involved are very small compared to the other terms. Equation (2.6) with the 
continuity equation for incompressible fluids which states Q = A,Vj = A 2 V 2 allows for the 
solution of this equation in terms of one of the velocities. Assuming that P ( > P 2 , i.e. the 
flow is from P, to P 2 , then solve for V, from the continuity equation as: 

H-4v, 

1 D , 2 

and substitute this velocity into Eq. (2.6) and solve for V 2 : 


V 2 =± 


2 (P l -P 2 ) 

P 

r, a*) 

l d;) 


(2.7) 


When the flow reverses, i.e. P, < P 2 , then the velocity at point 2 is described by the above 
equation with the pressure terms switched and a negative sign on the square root. The ideal 
volumetric flowrate (Q ideal ) for an incompressible fluid can be expressed as Q idea , = A*V. 

In a realistic flow situation though, there is a loss due to the Vena Contracta effect. This 
loss is empirically quantified by a discharge coefficient (C d ), which represents the 
percentage of the ideal flow that actually occurs. This coefficient, when multiplied by the 
ideal flow, yields Q^,, as: 

Q«„ = 0,0*, = AC„V (2.8) 

Substituting Eq. (2.7) into Eq. (2.8) for velocity: 


Qrcal ~ AC j 


< D, 

\ D lJ J 


■\[f{ ^ ~ * 4 > ^2 


(2.9) 


For the landing gear shown in Fig. 2-2, there are two flows that are of concern, the flow 
through the orifice plate and the flow into and out of the snubber chamber. Define Q s c as 
the flow rate into the snubber chamber in the compression mode, where the snubber orifice 
area (A s ) becomes A s c , which allows larger flow. The flow rate through the snubber 
orifice during the extension mode is defined as Q S E , and the area A s becomes A S E , which 
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only allows small, restricted flow. In both cases, the flow through the main orifice plate is 
Qo- 



Figure 2-5: Control volume between piston and orifice plate. 

Figure 2-5 shows the direction of fluid flow into and out of a control volume in the lower 
chamber as a function of stroke mode (extension or compression). In relating the flow 
rates to the pressures, defining a control volume as shown by the dashed line in Figure 2-5 
is necessary. The stroke rate is defined as: 

X s = X wg -X a (2.10) 

where the compression mode is given by X s > 0.0, and the extension mode by X s < 0.0. 
The flow is assumed to be negative leaving the control volume, and is positive entering it. 
For an incompressible fluid, the volumetric flow rates for compression and extension can 
be written as: 

Qo + Qs+^X, .= 0.0 ( 2 . 11 ) 

during the compression mode, and: 

Qo + Qs+A l X s = 0.0 (2.12) 

during the extension mode. Equation (2.9) defined the general form of the equation for a 
flow rate. Substituting the appropriate pressures, areas, and diameters into Eq. (2.9), the 
flow rate through the orifice plate during the compression mode (when flow is out of the 
control volume and P L > P u and P s ) can be written as: 
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Qo = -\c d 


(2.13) 
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for P L >P U 


where d 0 is the effective diameter of the main orifice, D L is the diameter of the lower 
chamber, and C d is the discharge coefficient of the main orifice. The flow through the 
snubber orifices during this mode is described by: 


4P^P S for P L >P S (2.14) 


with d s c as the diameter of a snubber orifice, D L as described above, is the discharge 
coefficient of the snubber orifice and A s c is the effective area of the snubber orifice. 
Similarly, for the extension mode, where flow is into the control volume (P L < P u and P s ): 




\S '-'dS 


( jC\ 




1- 


A. 


Qo = 4A 


2 


r, 

(d ] 

4N \ 

i p 

i— 



f 

l 

[d l j 

) 


~ f0r P u > P L 


(2.15) 


where the difference between this equation and Eq. (2.13) is that the pressure terms have 
exchanged positions and the whole term is now positive. The flow rate through the 
snubber orifices during the extension mode is given by: 


Qs = A E S C. 


dS 


' (d^ 

1- 


\ d rJ ) 


for P s > P L 


(2.16) 


where D R is the effective diameter of the annulus snubber chamber, d s E is the diameter of a 
snubber orifice, A S E is the effective area of the snubber orifices and C dS E is the discharge 
coefficient of the snubber orifices in the extension mode. To simplify Eqs. (2.13), (2.14), 
(2.15), and (2.16), let the non-pressure terms be redefined as: 
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respectively. Substituting Eqs. (2.13) and (2.14) into Eq. (2.11) and Eqs. (2.15) and 
(2.16) into Eq. (2.12) using this new notation, rewrite Eqs. (2.1 1) and (2.12) as: 


4K ~ p u ~ E i - p s + = 0.0 for X s > 0.0 (2. 1 la) 

4?u ~ P l + E 4 4P s ~ P l + = 0 ° f or X s < 0.0 (2. 12a) 



Figure 2-6: Control volume for the snubber chamber. 

Additional information about the flow rate-pressure relationship can be gained by 
studying a control volume in the snubber chamber as shown by the dashed line in Figure 
2-6. The variables A R and D R in Fig. 2-6 are the rebound chamber annulus area and 
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effective diameter respectively. P s is the pressure in the rebound chamber and d s c and d s E 
are the diameters of the snubber orifices in the compression mode and extension mode 
respectively. In the case of compression, where X s > 0.0 and P L > P s , 

<2f + A*X, =0.0 (2.17) 

Substituting the flow rate Qs C of Eq. (2.14) into Eq. (2.17) yields: 


-A c C c 

n s u dS 


f 


1- 


.A.J 


4K 


■p s +a r x s = 0.0 


(2.17) 


From previous notation of E,, this expression becomes: 


-E 2 ^^P s +A k X s =0.0 (2.18) 

Rearrange Eq. (2. 18) to get an expression for the pressures in terms of the stroke rate as: 


Ad 


t 2 

Substitute Eq. (2. 19) into Eq. (2. 11a) and solve for the variable P L as: 


(2.19) 




A, A r 
% 


\ 2 


V-2 


( 2 . 20 ) 


where P u is given in Eq. (2.3). Square both sides of Eq. (2.19) and solve for P s as: 


It)’*-' 


Ps = P\~ 

Similarly, for the extension case with X s < 0.0: 

P^Pu- 


( 2 . 21 ) 
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P, = Pl + 
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( 2 . 22 ) 

(2.23) 


These known pressures [Eqs. (2.3), (2.20), (2.21), (2.22), (2.23)] can now be substituted 
into Eqs. (2.1) and (2.2). Algebraic simplification of these equations leads to the 
compression and extension cases in terms of readily measurable quantities as: 
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M„K ~ M u g L + (A A l }P sj 
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(2.2a) 




for the compression case, and: 
MX, = M.g-L + {A,-A L )P„ 
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)A ! -W 


for the extension case. Introduce a new notation using subscripts to simplify the above 
equations: "1" and "2" will be associated with compression (equation set (a)), and "3" and 
"4" with extension (set (b)). With this change, the equations can be written in the form: 

MX S = M„g-L + C in X s 2 +K }n xr±f (2.1c) 

AAA = M L g+C 2l4 X s 2 + K 2/i Xr-F l +f (2.2c) 

where the coefficients of the stroke rate squared term are assigned the C-s, and the 
coefficients of the stroke position term are the K/s. 
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The only unknown term left in these equations is friction. As mentioned 
previously, friction in this gear comes mainly from two sources, friction due to tightness of 
the seal and friction due to the offset wheel (moment). The seal friction is assumed to be a 
maximum value statically and some function of velocity in the dynamic state. The 
functional relationship between frictional force level and velocity could be determined 
through testing. The friction due to the offset wheel is the result of the moment produced 
by the nonaxially loaded piston within the cylinder. 



Figure 2-7: Schematic of gear for friction model development. 

It can be seen from Figure 2-7 that the force between the piston head and the cylinder, N, is 
a result of the tire force, F t , applied at moment arm, ma, from the centerline of the piston. 
The frictional force due to the offset wheel (F ow ) is assumed to be of the form (refer to 
Figure 2-7): 

F 0W =pN (2.24) 

Where N is the normal force of the cylinder wall resisting the side of the piston head, and 
(1 is the coefficient of friction between the two parts. To find the unknown force N, sum 
the moments about point O to zero to get: 

ZM 0 : F t ma - N(X S + stp) = 0 (2.25) 

Where stp is the minimum distance between the piston head and the lower seal when the 
gear is fully extended. Rearrange Eq. (2.25) by isolating N, and then substitute N into Eq. 
(2.24) to get an explicit form of F ow : 
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N = 


ma* F t 


-C* -X m +stp 

F = Li 

aw r* 


ma* F, 


X w -X a +st Pj 


V 


The total friction in the landing gear, f, in equations (2. lc) and (2.2c) is now assumed to be: 

f-F- + F«, ( 226 ) 

This development assumes that a proportionate part of the fuselage (half of the 80% 
of the total weight that rests upon the main gear) is treated as a lump mass centered at the 
centerline of the main upper cylinder. Also, this model takes into account only vertical 
loads on the strut. The tire is modeled as a nonlinear spring and damper. This tire model 
does not take into account spinning stiffness (because the test tire does not spin) or spin-up 
drag. The fluid is assumed to be incompressible and all structural members are assumed 
to be rigid, with each having only a vertical degree of freedom. These assumptions are 
good only for straight-line taxiing over runway profiles and landing impact (spin-up drag 
on the tire does not significantly effect the vertical loads on the strut). Any braking or 
turning maneuvers are not covered in the development. The equations developed here are 
the basis for a "rollout" simulation. 


2.4 Summary 

In this chapter, the nonlinear equations of motion were developed for a general, 
telescoping main landing gear. These equations contain a pneumatic spring that is 
determined based on the polytropic gas compression law, a hydraulic damping that is 
proportional to the stroke rate squared, gravitational forces, lift, inputs from a runway, and 
finally friction, which is composed of both a constant seal friction and a variable bearing 
friction. These equations explicitly contain the empirical parameters of polytropic gas 
constant, discharge coefficients for both the main orifice and the snubber orifices, and the 
friction levels in the gear. These parameters are the only variables that appear in equations 
(2.1) and (2.2) that cannot be directly measured. 

Equations (2.1) and (2.2) are highly nonlinear and are discontinuous due to the 
differing values of friction and discharge coefficient as a function of extension and 
compression. Chapter 3 will discuss more about the nature of these equations and present 
a method of solving these equations for gear displacements and velocities. 
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Chapter 3: Numerical Analysis 


3.1 Introduction 

In Chapter 1, a brief history was given in regard to past and concurrent landing 
gear research. In Chapter 2, the terminology associated with a landing gear was defined 
and the equations of motion were developed. As seen in Chapter 2, the equations of 
motion of the landing gear system are nonlinear, due to velocity squared damping, 
polytropic spring rate, and friction. Many numerical routines exist to integrate nonlinear 
equations. However, those given by Eqs. (2.1) and (2.2) require some special 
consideration in that they are nonlinear, and discontinuous due mainly to friction, and 
under some conditions, stiff. This chapter will discuss the problem of stiff equations and 
discontinuities. It will detail the process undertaken to numerically integrate these types 
of equations and present a final scheme that successfully solves the problem. 

3.2 Model Integration 

The linearized model as presented by Edson and Ross 5 indicated that landing gear 
systems are stiff. In a system of two or more bodies, one type of "stiffness" in the 
equations is a result of the difference in time scales of the dynamics between the various 
bodies, or, in other words, there is at least one body whose solution time scale is much 
smaller or larger than the others 15 . The problem this causes for a numerical routine is that 
the time steps attempted must be small enough to accurately track the progress of the 
short time scale solution. These time steps can be orders of magnitude smaller than the 
time step required to accurately predict the solution of the other masses. The integration 
routine is therefore spending a lot of time tracking this fast solution while carrying other, 
slower solutions along. Most common integration routines are based on forward Euler 
schemes, known also as explicit integration routines. The problem of taking large time 
steps with an explicit Euler integration routine, when a fast time scale is present in the 
solution, is that the total solution may become unstable. However, a different type of 
routine, called an implicit (or backward) Euler routine, is able to take larger time steps 
when a fast time scale is present. These implicit routines are what are used to solve stiff 
equations. The main difference between the two is that the explicit routine uses derivative 
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information at the previous step to make the next step, whereas the implicit routine uses 
derivative information at the attempted step to reach that step. This assures that the 
implicit routine is stable, even for very large time steps, whereas the explicit routine is 
not. 

The equations of motion of the landing gear, as developed in Chapter 2 are 
numerically stiff. The numerical stiffness in the landing gear case comes from a couple 
of sources. Part of the stiffness is a consequence of the difference in time scales (about 
30 times difference) between the lower mass, which, because of its smaller mass 
compared to the upper mass, experiences very high accelerations, and the upper mass, 
which is very large and has much lower accelerations. The other source of numerical 
stiffness is introduced from the sliding friction model (Eq. 2.26). In the mathematical 
model, as velocity changes sign, the friction essentially steps from one large value to the 
negative of that value. In the original simulation, this was modeled as a discontinuous 
process, assigning a negative sign to friction as velocity passed through zero. To 
allieviate this discontinuity, a hyperbolic tangent was used as a continuous function to 
cause friction to change sign. It has been specified that the sliding friction will go from 
one value to the negative of that value in a velocity band around zero of about +/- 1 in/sec. 
This is a fix to the discontinuity problem, but introduces a new stiffness problem. The 
rise time of this function is taken to reduce numerical stiffness while maintaining a 
friction model that approaches reality. 

Many numerical routines were investigated in an attempt to solve this problem of 
stiff differential equations. It was found that a Runga-Kutta routine, with strict tolerances 
could solve the problem, but the time of solution was unacceptable. After further 
investigation, it was found that Eqs. (2.1) and (2.2), could be solved much more 
efficiently with an implicit predictor corrector routine. A predictor corrector routine uses 
a polynomial based on previous solution points to first extrapolate the solution to the next 
time step (predictor) and then uses correction iterations to drive the error between the 
predicted solution and the solution which satisfies the differential equation to within some 
tolerance 16 . It should be noticed that the predictor corrector routine needs some initial, 
one-step integration method to accumulate the first few points upon which to build the 
initial polynomial. For this landing gear case, a modified Adams-Moulton method was 
used. The routine, DDRIV2.f, was written at the Los Alamos National Laboratory (by 
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D. Kahaner and C. Sutherland, 9/24/85, available at http://gams.nist.gov/). The routine is 
designed to solve n first order ordinary differential equations in state space form given the 
initial conditions. The program also has options to allow the solution of both stiff and 
non-stiff differential equations, as well as an option to allow a dynamic selection of 
stiffness. For stiff equations, it uses a fifth order predictor-corrector and for non-stiff 
equations, it uses a twelfth order predictor-corrector. The routine is very flexible and 
contains checks to ensure proper usage. This program also has many input parameters 
that need to be selected with care, depending upon the problem to be solved. These 
parameters include the maximum time step attempted by the routine, defined by the 
difference between the initail time and the requested final time (0.00025 is used), a value 
of the requested relative accuracy in all solution components (le-6 for this problem), the 
smallest physically meaningful value for the solution (le-15), and the mode of stiffness 
solution (dynamic selection). The parameters in the current simulation have been set to 
values that seem to most efficiently solve the problem. 

As mentioned earlier, the problem at hand is also discontinuous. Two reasons 
exist for this discontinuity. The first occurs in the damping coefficients, Q's as presented 
in Eqs. (2.1c) and (2.2c). The damping coefficient is a function of fluid density, p, gear 
areas and discharge coefficients, C d 's, only. The discharge coefficients are assumed to be 
functions of orifice geometry (diameter) only. This model assumes that the flow through 
the orifice will be laminar (below a certain Reynolds number). A representative equation 
for the discharge equations is given 17 to be: 

C d = 0.8/? 2 -0.4813)3 + 0.8448 

In this model, p is the ratio of the orifice diameter over the diameter of the chamber from 
which the fluid is flowing. This model is for circular holes with rounded edges and was 
selected as a first approximation to the actual, unknown, discharge coefficient. If the gear 
is going from an extension to a compression mode, or vice versa, the value of the 
diameter of the rebound chamber inlets change nearly instantaneously. This is due to a 
slip ring in the physical gear that responds to the flow of the fluid and either chokes the 
flow (extension) or slides to a position that allows easier flow (compression). The model 
of this process is discontinuous. For compression, one value is used, and for extension, 
another value of discharge coefficient is used. Even though this is a discontinuity, it is a 
minor one. Since this discontinuity effects the calculation of a fluid damping coefficient 
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that gets multiplied by a velocity squared (see Eqs. (2.1c) and (2.2c)) term, and this 
discontinuity occurs at zero velocity, the effect of this discontinuity on the solution is 
small, and no further steps were taken to smooth the transition. 

The second discontinuity comes from the event of strut sticking or breaking loose. 
The event of stiction, or the sticking together of the two parts to move as a rigid body, is 
modeled in this simulation. The method for implementing stiction friction, as used in the 
simulation, was developed by Kamopp 18 . This model treats near zero relative velocity 
stick friction in a manner that does not introduce further numerical stiffness and does not 
require reformulation of the equations of motion. 

3.3 Kamopp Friction Model’ 8 

This section deals with the manner in which a numerical integrator can handle the 
task of integrating a model that includes a frictional model that allows sticking. Figure 3- 
1 shows a two mass system in which W,, Vj, M,, F, are the momentum, velocity, mass 
and applied force to the first mass. The same holds for the second mass. F r and V r are 
the relative force and velocity between the two bodies. 



Figure 3-1: Simple two mass system with stick-slip friction. 


For this two mass system, take the momentum (W ; ) of each mass as the state vectors. 
The state equations are then given as: 

Wt=F t -F r (3.1) 

W 2 =F 2 + F r (3.2) 

and the velocities can be solved from the momentum to be: 


26 





(3.3) 


Vm R 

1 M, 

K-i 


with V r = V, - V 2 . When the two masses are stuck together there should be no relative 
velocity. Let F r = F slick , the force required to keep V r = 0. For the two masses to display 
no relative motion through time, the time derivative of the relative velocity also needs to 


be zero. This derivative is taken as: 


y = R_R 

r Mj M 2 


Substitute Eqs. (3.1) and (3.2) into Eq. (3.5) to get: 

(Jj-Jy) (5+e) 


Setting the Eq. (3.6) to zero and solving the resulting expression for F r , or F stick , gives: 

F. = — — f; ^ — F 2 (3.7; 

s " ck m,+m 2 m,+m 2 

The logic of Kamopp's model states that if the absolute value of the relative velocity, IV r l, 
is smaller than some defined quantity, 5, and if the absolute value of the difference 
between applied forces, IF,-F 2 I, is less than the peak, sticking frictional force, F^, then 
the masses will stick together. Otherwise, slippage occurs and there is relative velocity 
between the two masses, and friction can be any arbitrary function. For the case of a 
landing gear, the equations of motion can be cast into a form in which momentum is the 
state vector as: 

= MX, = M u g-L + CX-K U3 X: r -Fr (2.i: 

W 2 = M,X a = M,g + C 2 X] + K 1/4 x; r -F' + Fr ( 2 . 2 ] 

where, F r = +/- f. To use Kamopp's model, let: 

F\ = M u g - L + CX + K\nX~ r 

This reassignment allows the Equations (2.1) and (2.2) to be written in the form: 

W 2 = F 2 +F r ’ 
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after which the same logic as above can be applied. 

This model uses the relative velocity and the relative force as the decision factors 
for sticking. If the velocity is very close to zero Md the relative force is below a 
preselected sticking frictional force, then the logic is to assign the frictional force such that 
the relative acceleration is zero. This implies that the relative velocity remains constant, at 
some small value, 8, which is unwanted. An addition to this model was to damp out this 
small velocity. Coulomb damping was used to decrease the remaining velocity to zero 
after the sticking condition is active. When slipping, the friction can be any arbitrary 
function. 

In considering how to model friction, two other models were also investigated 19 . 
One was called the "bristle model" in which a number of bristles, or springs, are defined 
between the two relative surfaces. Each bristle has a stiffness and can be broken after a 
certain amount of relative force has built up. As some bristles break, others are 
established. The number of bristles established is a function of velocity. This model will 
capture the effect of sticking, but is numerically inefficient. The other model is called the 
"reset integrator model". This model is similar to the bristle model except that there is 
only a single bond between the relative surfaces. Its advantage is that it also has a 
mechanism for the frictional energy to damp out when the two bodies are sticking. This 
model is much more efficient than the bristle model and compares well to the Kamopp 
model. The Kamopp model was chosen because of its ease of implementation, its 
realism in capturing the slip-stick phenomenon, and its numerical efficiency. 

3.4 Treatment of Discontinuities 

For the reasons due mainly to stick-slip friction, the simulation contains 
discontinuities that need special treatment. The Adams-Moulton integration routine 
incorporates a system of warnings and errors that allows the user to become aware of 
some of the problems that the routine is having. One such warning to the user indicates 
when DDRTV2 is attempting too many iterations, or reductions of time step, to get to the 
next output time with the specified accuracy. This warning has been used to determine 
when a discontinuity, either sticking or breaking-free, has been encountered. The 
predictor-corrector routine is trying to fit the next point of at least a fifth order polynomial 
to a comer in the solution history (the discontinuity). The calling program has been 
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modified so that when this warning is activated, the main program switches to an error 
control, variable-step fourth order Runga-Kutta integration routine, RKF4.f (by S. 
Baudendistel and G. Haigler, 4/1/83, available at http://gams.nist.gov/). This routine is an 
explicit-type one step integrator and is based on Fehlberg's formulas. This program was 
used to get past the discontinuity, at which point, the main program directs the predictor 
corrector to continue the solution. The variable step feature of this R-K routine is useful 
in error control. When passing the solutions from Adams-Moulton to R-K, and back 
again, it was found that the error tolerance of each program needs to be near the same 
order. Numerical errors in the form of instabilities were encountered when different 
tolerances were used. The solution is unstable for a difference in specified tolerance of 
three orders of magnitude, i.e. R-K tol. = le-3 and A/M tol. = le-6. For error tolerances 
within two orders of magnitude, the over all solution was stable, but there were many 
areas of local numerical instability. Finally, for R-K tolerances of le-5 and A/M 
tolerances of le-6, the solution is well behaved, with only a few numerical problems 
under stick-slip conditions. It was found that decreasing the R-K or AM tolerance to le- 
7 was a bad trade off between the time it takes to complete the run and the incremental 
increase in numerical stability. When the gear breaks loose, the predictor corrector will 
generally call Runga-Kutta once or twice, depending on how quickly the break-loose is. 
The faster the break loose, the less it calls R-K. Adams-Moulton seems to have the most 
problem when the gear goes from a relative motion state to the stuck state. This condition 
will almost always trigger the R-K. When the gear experiences forces and velocities very 
near the break-free point, i.e., it is in a continuous state of sticking and slipping, the run 
times can become longer, as Rung-Kutta does more of the integration. Under the current 
set of parameters, the run times for a fully dynamic case (no sticking) is about 25 seconds 
real time per 1 second simulated time. When sticking and slipping are involved, the run 
times are longer, about 2 minutes real time per 1 second simulated time. 

3.5 Summary 

In an effort to maintain model fidelity, the equations were left in their nonlinear 
form rather than linearized. Such considerations as the stiffening effect of sliding friction 
and the discontinuous behavior of stick-slip friction and discharge coefficients were also 
factors in this decision. Therefore, numerical routines were found to handle this problem. 


29 



A modified Adams-Moulton routine integrates the stiff, nonlinear equations until a 
discontinuity is encountered, as detected by the routine attempting to reduce the time step 
too many times without getting to the next step, at which point, a variable step, error 
control Runga-Kutta integrates past the discontinuity, and then the Adams-Moulton 
continues the solution. Adjustment of the control parameters to the integration routine 
has led to a more stable solution as well as reasonable run times. The next task is to 
verify the model parameters with experimental data. Chapter 4 will detail the facility and 
equipment used in the tests and Chapter 5 will present experimental results and discuss 
their significance and usefulness in the validation process of the simulation, as well as 
present dynamic comparisons between the updated model and test data. 
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Chapter 4: Experimental Facility 


4.1 Introduction 

The equations of motion as developed in Chapter 2 are nonlinear, due to the 
velocity squared damping term and the polytropic gas law assumption, and stiff and 
discontinuous due to friction. As discussed in Chapter 3, however, a number of 
numerical integration schemes were evaluated for use in this simulation. The final method 
uses a predictor corrector and a Runga-Kutta scheme to solve the problem. This chapter 
describes the experimental facility and equipment used to validate the simulation with 
experimental data. 

The objective of the testing was to determine the physical characteristics of the 
A-6 gear and to use that information to adjust parameters and/or models in the 
simulation. Quasi-static tests determined such quantities as masses, maximum static 
frictional forces, load-stroke curve for the nitrogen spring, and tire load-deflection curve. 
Dynamic tests were used to find dynamic levels of friction and values of orifice discharge 
coefficients. Initial tests to validate the simulation software were performed at NASA 
Langley Research Center. The particular equipment used was an instrumented A-6 main 
landing gear and a mobile data acquisition system. The gear is mounted on a truss-like 
drop carriage, which is constrained for vertical motion within a main, translational 
carriage. The tire of the gear rests on a hydraulic shaker table which is controllable via 
computer. 

4.2 Test Equipment 

As stated previously, an A-6 main landing gear was selected for these tests. This 
gear was chosen for its availability. It and four other main gears were scrapped by the 
Navy as part of the phasing out of the A-6 Intruder fleet. The landing gears are still in 
operational condition and were acquired from NAVICP-PHILA, a Naval surplus yard, as 
a gift toward research. The gear and a GoodYear USA 36X1 1 Type VII tire inflated to 
120 psi was installed on the drop carriage so that it would be in the standard vertical 
position, as shown if Figure 4-1. A connecting plate was fabricated to allow the normal 
mounting of the gear to the plate, and the plate was then rigidly connected to the drop 
carriage. The drop carriage is a truss-structure that weighs about 4.5 tons and allows the 
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gear to be raised and lowered. The translational carriage weighs about 55 tons and rides 
on horizontal tracks. It can be moved such that the landing gear tire is over concrete only 
(for drop tests) or over the shaker table (for some static and many dynamic tests). The 
mass of the drop carriage rests upon the landing gear. This mass simulates the rigid 
portion of the aircraft mass carried by the gear. Once the gear is loaded, the shaker table 
is used to input forces into the gear. Hydraulic lift cylinders, powered by a hydraulic 
mule, are used to lift the drop carriage and unload the gear. Once the gear has been lifted, 
the ability exists to lock the gear in that position with hydraulic valves. 



Figure 4-1: Schematic of experimental set-up. 


The hydraulic shaker table was built specifically for the task of examining the A-6 
landing gear. It was built by TEAM Corporation to the specifications of NASA LaRC. 
These specifications included the capability to perform a step bump of one inch in no 
longer than 2 ms while bearing 12,000 lbm. The shaker is also capable of simulating 


32 





wave functions at user-selected frequencies with amplitudes of about 3.5 inches, at a 
dynamic force level of at least 10,000 lbf., and for a duration of at least two cycles within 
no more than a ten second period. The wave functions include: (1 -cos), sine, a 
trapezoidal bump with user-selected rise time, and a saw-tooth wave form. The shaker 
can also be driven by a file containing runway elevation versus time data and, through 
positional feedback to the controller, internally adjust the inputs to the shaker to 
accomplish the input profile. The shaker is also capable of supporting variable static 
loads of at least 12,000 lbf and allows actuator movement of 6 inches. This shaker 
package included a digital servo control system that operates from a PC computer. This 
controller provides for user-selectable displacement, velocity, or acceleration actuation of 
the shaker head. It is also capable of controlling the shaker to accomplish all of the built- 
in waveforms and user selected profiles. This software also provides plots to show the 
user-selected runway profiles/simulations versus the accomplished runway 
profile/simulation. 

The gear was instrumented to provide the necessary information for model 
validation (see Fig. 4-2). There are two accelerometers, one placed at the upper mass and 
the second one at the lower mass. Two potentiometers are also used, one to locate the 
upper mass with respect to a fixed position on the translational carriage and one to 
measure the relative position between the upper and lower masses of the gear. Two 
pressure transducers are included in the instrumentation as a check of some of the basic 
assumptions of the simulation (mainly that the fluid and the gas do not mix to any 
significant degree after initial shaking). One is located just outside the charge port of the 
upper cylinder, and the other is embedded in the piston head. Finally, there is a strain 
gage on the wheel axle of the gear. This gage is calibrated to read the vertical load 
through the strut and the bending moments induced by the tire. 

Table 4-1 shows the instrument sensitivity and other detailed sensory information. 
These instruments were selected to allow direct comparisons to the simulation results 
developed from the equations of motion obtained in Chapter 2. 
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Figure 4-2: Instrumented A-6 landing gear. 


NAME 

Type 

Range 

Offset 

Sensitivity 

*** 

*** 

*** 

Co 

Cl 

Upper Mass Position 

Slide Pot Wire, TCC 

40 in. 

28 in 

-10.48 in/volt 

Strut Piston Position 

Slide Pot, Bourns 

16 in. 

-1.94 in 

4.04 in/volt 

Lower Chamber Press. 

Pressure Transducer, Kulite 

2 ksi. 


Upper Chamber Press. 

Pressure Transducer, Kulite 

2 ksi. 


Upper Mass Accel. 

Accelerometer, Kistler 

Kama 

0 


Lower Mass Accel. 

Accelerometer, Kistler 

itaitna 

0 

11^23325181 


Wire Strain Gage, MMT 

162.2 klbs. 

1050 lb 

6488.0 lb/mvolt 


BLH 20 klbs 

20 klbs. 

161b 

3999.20 lb/volt 

| Shaker Head Position 

LVDT, TEAM 

(+/-) 3.89 in. 

0 

0.77 in/volt 

Engineering Units (EU) 

= Co + Cl * Voltage Reading 





Table 4-1: Instrument guide on A-6 test specimen. 
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A mobile data acquisition system has been developed to gather, manipulate, plot 
and store data taken from the tests. This system is a roll-around rack that allows 16 
channels (expandable) of input and incorporates a LAB VIEW interface. Data from two 
channels can be plotted in real time. At post test, up to 16 channels can be plotted versus 
time simultaneously, or any selected channel can be plotted against any other channel. 

The system has a user defined acquisition rate of between 1 Hz and 3000 Hz and has 
built-in user selected digital data filters. Finally, this system allows manipulation of the 
data and will store the data in a Microsoft EXCEL worksheet format. 

4.3 Summary 

The object of the tests, again, is to determine the physical characteristics of the A- 
6 test gear and use that information to update the simulation. These tests to validate the 
simulation software were performed at the Aircraft Landing Dynamics Facility in 
building 1262 at NASA Langley Research Center. An instrumented A-6 main landing 
gear is mounted on a truss-like drop carriage, which is constrained to vertical motion 
within a main, translational carriage. The tire of the gear rests on a hydraulic shaker table 
which is controllable via computer. A mobile data acquisition system records and 
manipulates the data incoming from the test set-up. Chapter 5 will explain the procedures 
of each test and present the results, as well as explain how the results of each test are to be 
incorporated into the model. 
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Chapter 5: A-6 Experimental Parameter Determination 


5.1 Introduction 

The previous chapters have defined the theoretical basis for modeling the landing 
gear, the numerical analysis involved in solving the equations of motion, and finally the 
test equipment used to determine some of the unknown parameters of the model. This 
chapter is divided into two sections. The first section describes the procedures and results 
of the tests to determine static values of some parameters like system masses, static 
frictional levels, pressure-stroke curve, and tire load-deflection curve. The second section 
describes the procedures for determining dynamic parameters in terms of discharge 
coefficients and dynamic frictional levels, and tire damping. Adjustments to the statically 
updated model are made using dynamic data and the final model is compared in 
frequency space to test data. 

5.2 Determination of Static Parameters 

The first set of quasi-static tests were designed to define the masses of the system 
and some static frictional loads. It was recognized early on that there may be external 
forces acting on the upper mass due to the friction in the bearings of guide rollers on the 
drop carriage. The first test was designed to measure the total system mass and to 
measure the frictional level of the bearings. This test was a quasi-static test and was 
performed as follows. The upper and lower mass were locked together to prevent relative 
motion and the drop carriage was raised by the lift cylinders until the tire lost contact with 
the shaker head. A load cell was then placed under the jack lug of the gear, shown in 
Figure 5-1, and the drop carriage was lowered until the entire mass rested on the load cell. 
Very slowly, the shaker head was manually raised and lowered for just over one cycle for 
a total displacement of about six inches. During this time, load cell reading and upper 
mass position were being recorded. The expected result was a hysteresis loop centered 
around the weight of the system with the loop defining the positive and negative range of 
maximum sticking friction in the carriage bearings. The test results are shown in Figure 
5-2. The horizontal lines describing the means of friction and weight on the plot were 
found by first, summing the entire load array and dividing by the number of points, thus 
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getting an estimate of value of the center of the loop (some of the data was repeated as the 
return stroke overlapped previous distance traveled and so the average is weighted toward 
the lower bound). Then, the data set was divided along this value into those values higher 
than the average and those lower. A mean was then found for each of these data arrays, 
defining the frictional upper and lower limits. These two values were then summed and 
the average taken to find the average total system weight to be 9465 lbs. Compared to the 
total weight of the system, the static frictional level of the guide rollers of 1 17.7 lb. is only 
1.25% of the load felt through the gear. Under dynamic conditions, this frictional level 
will decrease, having an even smaller effect on the dynamics. For this reason, this 
external force is neglected in the simulation. It's addition would complicate the Kamopp 
model of friction and add only a very small increase of fidelity. 



Figure 5-1: Load cell under jack lug to measure system mass and friction. 
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Weight of System and Friction in Carriage 



Figure 5-2: Total weight of the system and frictional hysteresis loop. 

The next test dealt with finding the mass of the piston, the wheel and tire, and the 
fluid inside the piston. For this test, the upper chamber was vented to the atmosphere so 
that the air spring was taken from consideration. The load cell was under the jack lug as 
in the first test and the gear started from the fully compressed position. The lift cylinders 
were used to slowly raise the upper mass through about twelve inches of the gear's stroke 
and returned it to the fully compressed state very slowly. The expected result of this test 
was another hysteresis loop as found in the first test. However, the center of this loop 
would be the lower mass weight and the boundaries would be the constant seal static 
friction. Figure 5-3 displays the test results as measured in the lab, showing the weight 
to be 318.4 lbs. and a frictional band of +/-1 15.7 lbs. The lines indicating the mean 
weight and the mean values of friction were obtained in the same way as in the first test. 
A check of the test accuracy was also performed. The wheel and piston of one of the 
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extra gears and a gallon and a half of fluid were measured on a scale. The total lower 
mass as measured by the scale was 320 lbs. This agrees very well with the data found in 
the quasi-static test. 


Weight of Lower Mass and Seal Friction in Strut 



Figure 5-3: Weight of lower mass and frictional hysteresis loop. 

The third experiment was also a quasi-static test. The goal of this test was to 
obtain data concerning the tire load-deflection relationship. The gear was serviced by 
checking the fluid amount and adding nitrogen to the upper chamber to get a static stroke 
of about 2 inches above what is desired. The strut was then exercised through vigorous 
shaking via input from the shaker table. A steady state was reached between the fluid and 
the compressed gas, i.e. some of the gas volume was lost due to being dissolved into the 
fluid, and the gear settled to near the desired static stroke. A procedure was developed 
along these lines to try to predict what charge pressure to inject to a fully extended gear 
that has not yet been shaken. It was found that if the goal of static stroke was 3.5 inches. 


39 




for example, friction could cause the gear to stick at a static stroke value as high as 6.5 
inches. With a little shaking, the gear could be settled again to 3.5 inches, indicating that 
friction will significantly affect procedures in the quasi-static regime. The test started with 
the gear fully extended and the tire above a platen which rested on the load cell. The lift 
cylinders were used to slowly lower the gear until it came to equilibrium and then to raise 
the gear again. As a means to get more data, the drop carriage was locked in position to 
allow no upper mass displacement and the shaker head was used to further deflect the 
tire. Two points of data were taken from this test. After some initial deflection, aircraft 
tire spring behavior becomes essentially linear. The data set found from the continuous 
test was combined with these two points. A third order polynomial was used to represent 
the data to capture both the nonlinear behavior at initial compression and the linear 
behavior around the operating point at 1.6 inches. Figure 5-4 shows the data and the 
cubic fit. 
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It may be noticed in Figure 5-4 that the tire load data does not start at zero when tire 
deflection is zero. The positional data of tire deflection from the continuous data was 
calculated from measured quantities (Tire deflection = X wg - X s ). The intercept of the data 
(using an average through the hysteresis) is approximately -0.17 inches. It is thought that 
combined calibration error of the upper mass (X wg ) and strut positional (X s ) 
measurements, as well as some flexure in the drop carriage may account for this 
discrepancy. 

A final quasi-static test was performed to gather information concerning the 
pressure-stroke curve and to help correlate the offset wheel friction model. For this test, 
the upper mass was locked into its equilibrium position by using tie down cables to 
prevent motion in the upward direction while the lift cylinders prevented motion in the 
downward direction. The shaker head was then moved from the zero point to the fully 
retracted position, allowing the gear to stroke about 7.25 inches. Two runs of this test 
were made. The first controlled the shaker head to move vertically at a slow rate of about 
0.084 in/sec till no further stroke was possible with the shaker head (stroke of about 2.8 
in left on the gear). The second test used a stroke rate of about 0.725 in/sec. These two 
tests were performed to determine the effect of nitrogen gas dissolving into the hydraulic 
fluid and to determine some of the effect velocity has on frictional levels. It was found 
(see Figure 5-5) that for a faster compression rate, less gas is dissolved into the fluid, 
leaving more gas in the chamber, causing the air spring to be stiffer. The hysteresis in the 
pressure measurements represents the amount of volume of gas lost to or gained from 
the fluid. No plans are made to model this effect. However, the decision as to which 
pressure curve to use as the model is made with the reasoning that stroke rates that are 
results of runway inputs will be high. Therefore, the air curve found during the higher- 
rate test was selected to represent the dynamic response of pressure to stroke. A curve 
using the form of Eqn. (2.3) was used to fit this data. This calculated curve is also shown 
on Figure 5-5. Two points in the extrapolated area of this curve were checked against test 
data. The first point, at about 1 1.0 inches, agreed to within 1.6% of the test data and the 
second point, at fully extended stroke, 15.09 inches, agreed to within 9.6%. 
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Pressure-Stroke Curve for A-6 Landing Gear 



Figure 5-5: Pressure-stroke curve and fitted analytical expression. 

The folowing procedures were developed in an attempt to statically quantify 
frictional effects. The result of this method was unclear and the final frictional model was 
developed using dynamic test data. The approach to finding the analytical expression for 
friction using static data is to use both the load data measured from the axle strain gage 
and to calculate the "theoretical" load (or frictionless load) by using the nitrogen pressure 
times the area it acts on. By subtracting these two data sets, the remaining loads are 
friction induced. However, as described in the previous section, the pressure does not 
follow a set rule because of the solubility of the gas into the fluid. It is therefore assumed 
that the frictional loads should be symmetric with the zero load axis. Figure 5-6 shows 
the result of subtracting the theoretical "pressure force" from the measured axle load. In 
this plot, two data sets are represented. The first represents the load measurements that 
were taken through a slow compression and extension around the 2.2 to 7.5 inch range. 
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The second represents data that were taken during a faster compression and extension rate 
through a stroke range from about 5.0 inches to 15.0 inches. There are two major notes 
to make from this plot. The first is the pressure effect as mentioned above. It is believed 
that from one time to the next, in a quasi-static regime, the pressure cannot be accurately 
predicted because of the solubility effect. So, even though the area is constant, the 
pressure for a given stroke value is very transient unless a long period of time is allowed 
for the gas and fluid to come to equilibrium, or the process is done so rapidly as to allow 
no mixing. This transience of pressure can explain why the two data sets in Figure 5-6 
are not centered around zero. The next point to notice from the plot is the difference of 
scale between the two data sets. The set taken at a much slower rate shows a much 
greater frictional hysteresis loop, whereas in the other set, where the stroke rate was 
faster, the frictional loop is thinner. This lends credibility to the theory that sliding friction 
is a function of velocity. Statically, one encounters the maximum amount of friction 
possible. As velocity decreases, the friction also decreases. Beyond some velocity value, 
the friction remains essentially constant. 


Pressure Corrected Axle Load-Stroke Friction Data 



Figure 5-6: Result of "pressure load" subtracted from axle load measurements. 
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In light of these two points, and the assumption that friction be centered (or 
symmetric) around zero, a process was developed to center the test data. A rough median 
value of the two data sets in Figure 5-6 was found and those values were subtracted from 
their respective data sets, resulting in Figure 5-7. 


Corrected Axle Load-Stroke Friction Data 



Figure 5-7: Zero centered frictional load data. 

The data, as seen in Figure 5-7, may seem to indicate a functional relationship with strut 
stroke and velocity, but to draw any conclusions from this data, in light of how the 
pressures are capable of change, would not be useful. In light of the uncertainty 
associated with this quasistatic frictional information, a working model of friction was 
developed using dynamic information only. The initial model for comparison against 
dynamic data will contain the frictional model as developed in Chapter 2, and static 
friction values of (arbitrarily) 1.3 times that of sliding friction. 
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In summary, the system weight was found to be about 9465 lbs. and the lower 
mass about 318 lbs., leaving the upper mass to be about 9147 lbs. A frictional level of 
+/- 117 lbs. was found to exist in the guide rollers of the drop carriage, but in comparison 
to 9147 lbs., this force is negligible. A static seal friction of roughly 1 15 lbs. was found 
to exist in the strut. The friction due to the offset wheel was found possibly to be a 
function of stroke and gear force, but the data is not clear enough to draw a firm 
conclusion. The pressure-stroke curve was found to best characterize a rapidly 
encountered stroke. Finally, the tire load-deflection curve was found to be represented 
best by a cubic equation to capture both the nonlinear effects at initial compression and the 
linear spring rate at normal operating deflections. These values and models, except those 
noted, were implemented for comparison to test data in a dynamic regime. 

5.3 Dynamic Testing 

Many tests need to be performed to fully validate this landing gear model, such as 
step bumps, ramp inputs, varying sinusoidal inputs, etc. The remaining unknowns, after 
the static testing, are the tire damping coefficient, a dynamic polytropic gas constant, the 
three discharge coefficients, levels of sliding friction, and criteria for when the strut breaks 
free from the static friction and starts to slide, and when the strut sticks again. A 
frequency response comparison between the test gear and the simulation of the gear to a 
sinusoidal sweep from a runway input was used to determine these quantities. This 
process is only to compare the gain amplitudes and phase shifts of varius parameters 
given certain, specified inputs to this nonlinear system. This process should not be 
confused with the typical frequency response used when dealing with linear systems. 

The purpose of these tests was to make a comparison between the model, which has been 
updated by the static data, and the test gear. These frequency response tests should 
demonstrate precisely how the landing gear system responds to a frequency sweep of 
given amplitudes. The variables under consideration, in comparing the two models, are 
the gear positional variables and pressure values in both the upper and lower chambers. 
The maximum, break away friction can also be observed by slowly increasing the 
frequency of a given amplitude and noting the force levels present in the gear when the 
strut starts to stroke. These variables considered in the frequency responses comparisons 
are important to verify because of the intended use of the program. The simulation will 
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be used to evaluate analytically the effect of various active control concepts that will be 
applied to the test set-up. 

For this frequency sweep test, the gear starts at rest on the shaker head. It is then 
shaken to allow the fluid and gas to come to equilibrium. The test used as an input a 
swept sine wave from 0.75 to 3.75 hertz in the course of 40 seconds with an amplitude of 
about 1 .0 inch. The model parameters of sticking friction, sliding friction, tire damping, 
discharge coefficients, and polytropic gas constant were adjusted such that the predicted 
frequency response and the test data were within about 10% agreement in the positioned 
and pressure variables over the whole frequency range. 

The method for setting the sticking friction was to observe the force level (above 
the weight of the system) at which the strut started to stroke. This quantity is about 2700 
lbs. for the 1 .0 inch amplitude case. The predicted slip time, as calculated by the 
simulation, is very sensitive to this number. The criteria for the break-away friction level 
may also be function of something other than purely force. From comparison with other 
runs, such as step bumps and other amplitude sine waves, this number is not consistent 
and may be a function of how long the strut has been stuck and how fast it was traveling 
when it stuck, and potential and kinetic energies in the landing gear system. No attempt 
was made to model these effects, and the constant value of 2690 lbs. was used in the 
following results. 

The tire damping in the model was adjusted such that the tire deflection predicted 
by the simulation and the data recorded in the test before the strut started stroking 
matched. The higher frequency tire mode was inspected in this region where only the tire 
is bouncing and the damping was selected to match the phase of the tire as well as 
possible. 

The other parameters of sliding friction, discharge coefficients and polytropic gas 
constant are coupled in that they all effect the stroke. The sliding friction and discharge 
coefficients work to damp the stroke and the polytropic gas constant affects the stroke by 
changing the curve of pressures as a function of stroke, making the spring stiffer or 
softer. These parameters were selected by inspecting the dynamic behavior of the 
variables of upper mass position (X wg ), strut stroke position (X s ), upper chamber, or 
pneumatic pressure (P u ), and lower chamber, or hydraulic pressure (P L ). If, for example, 
the predicted pressures were much larger than the measured pressures but the positional 


46 



variables seemed close, the polytropic gas constant needed to be decreased. However, if 
both the pressures and the positional data were greater than the measured data, then the 
damping in the strut, alone, may need to be increased. These parameters were adjusted in 
an iterative fashion until the frequency responses of each of the above mentioned 
predicted variables were within about 10% of the measured data. The value used for 
sliding friction is a constant 400 lbs. This value is somewhat arbitrary and in reality, may 
indeed be a function of stroke velocity. Further tests need to be run to quantify friction to 
a larger extent. 

The discharge coefficients are still fuctions of orifice geometry, as described in 
Chapter 3. However, a percentage of each of the three discharge coefficients is taken to 
help match the data. The model of the discharge coefficients includes only geometry 
information, no flow velocity information. These models may not be entirely appropriate 
for this landing gear, but the adjustments as mentioned above improve the results and are 
taken to be sufficient control on these coefficients. As long as the coefficients do not 
exceed one, and are greater than about 0.8, arbitrarily, these values are reasonable. The 
average of the discharge coefficients is about 0.9 

The polytropic gas constant was adjusted down to 1.1 from 1.19 as indicated 
earlier in this chapter. This change reduced the pressure amplitudes for given values of 
stroke. This change, along with the others mentioned previously, was sufficient to bring 
the predicted results to within 10% of the measured quantities. 

The frequency response comparison of both the measured data and the nonlinear 
simulation data are presented in Figures 5-8 (a-d). Figure 5-8a shows the gain and phase 
plots for the strut stroke variable (X s ) as a response to the 1 inch shaker head swept sine 
input. As can be seen from the figure, the maximum gain predicted by the simulation at 
natural frequency of about 1.6 Hz. is not as great as that recorded in the experiment by 
about 13%. This value is greatly affected by the damping in the strut. The parameters of 
discharge coefficient and sliding friction may still need to be adjusted. The phase of the 
response agrees very well with that of the test data. Over the rest of the frequency range, 
the two gains match to within an average of 5%. The "jumps" at the beginning and end 
of some of the phase plots presented are the usual difficulty in distinguishing between +/- 
1 80 degrees and the fast Fourier transformation routine calculating a response to a 
nonlinear system where no data exists. The next figure. Figure 5-8b, shows the 
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frequency response comparison of the upper mass position (X wg ) with respect to the 
shaker input. The gain plot shown here indicates that for higher frequencies, the upper 
mass, as simulated, moves through larger amplitudes (about 14%) than the test gear. 

One reason for this may be the frictional and discharge coefficients as mention above. If 
the strut were allowed to stroke a little further each time, the upper mass amplitudes 
would not be as large. The suggestion is to reduce the damping, to some small extent, in 
the strut to accommodate this behavior. The phase of this predicted value is also slightly 
in error. This may also improve with the suggested changes. 







Response of Wing/Gear Position due to Shaker Head (1 .0 in.) 




Figure 5-8b: Frequency response comparison of upper mass position to shaker input. 

Figure 5-8c shows the frequency response comparison of the pressure in the upper 
chamber (P u ). For the pressure measurements, it is useful to compare amplitude gains 
between the model and the actual measurements over a range of frequency inputs. The 
upper chamber and lower chamber (see Figure 5-8d) pressures are both high by about 
10% through all frequency ranges except near 1.6 Hz, the natural frequency of the gear. 
Again, damping in the strut is the critical component in this comparison. It is noticed that 
at natural frequency, where the strut would otherwise stroke more, the calculated pressure 
is not larger than the measured pressure. This lower pressure value corresponds directly 
to lower stroke values at this frequency. The damping limiting the stroke indicates that 
the main orifice discharge coefficient may need to be a little larger, to allow more flow 
through the orifice. This would allow larger stroke rate by decreasing the damping in the 
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strut. The polytropic gas constant may also need to be decreased slightly to lower the 
pressure amplitudes over the whole range of frequencies. 


Response of Pneumatic Pressure due to Shaker Head (1.0 in.) 



Figure 5-8c: Frequency response comparison of upper chamber pressure to shaker input. 


Response of Hydraulic Pressure due to Shaker Head (1.0 in.) 



Figure 5-8d: Frequency response comparison of lower chamber pressure to shaker input 
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Figure 5-8d of the hydraulic pressure shows a response similar to that of the pneumatic 
pressure. The difference in amplitudes (pressure) between the hydraulic and pneumatic 
drive the motion of the strut. 

5.4 Validation of Updated Model 

The swept sine used in Section 5.3 is an extreme case for a landing gear to 
encounter. The model parameters were selected such that the model predicted certain 
measured quantities to within 10% for this one case. As a validation to the model 
updates, several other tests were performed and the results were compared to the 
predicted data. The first was a swept sine test over the same frequency range and with an 
amplitude of about 0.5 inches, again over the course of 40 seconds. The results are 
shown in Figures 5-9 (a-d). The reasons, as discussed in the previous section, to reduce 
strut damping are apparent when looking at Figure 5-9a. Over the whole range of 
frequencies, the response as predicted by the simulation is below that of the test data by 
about 9%. With the exception of this one variable, however, the data and the predicted 
values agree very well. 


Response of Strut Stroke Position due to Shaker Head (0.5 in.) 




Figure 5-9a: Frequency response comparison of strut stroke to shaker input. 
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Response of Wing/Gear Position due to Shaker Head (0.5 in.) 




Figure 5-9b: Frequency response comparison of upper mass position to shaker input. 


Response of Pneumatic Pressure due to Shaker Head (0.5 in.) 



Figure 5-9c: Frequency response comparison of upper chamber pressure to shaker input. 
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Response of Hydraulic Pressure due to Shaker Head (0.5 in.) 



Figure 5-8d: Frequency response comparison of lower chamber pressure to shaker input 


The second check was to run a test where the sweep rate of the frequency range 
was different. This test swept a 1.0 inch amplitude sine wave from 0.75 to 3.75 Hz. over 
the course of 25 seconds. 


Response of Strut Stroke Position due to Shaker Head (1.0 in.) 









Phase (Deg.) 


4 


Response of Wing/Gear Position due to Shaker Head (1.0 in.) 



1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 

Frequency (Hz) 

Figure 5-10b: Frequency response comparison of upper mass position to shaker input. 



Figure 5-10c: Comparison of responses of upper chamber pressure to shaker input. 
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Response of Hydraulic Pressure due to Shaker Head (1 .0 in.) 



Figure 5-10c: Comparison of responses of lower chamber pressure to shaker input. 

As can be seen from Figures 5-9 and 5- 1 0, the simulation predicts the response of the 
measured values exceptionaly well (within 5%). The results to these past two tests, 
where the amplitude and then the sweep rate were changed, indicate that the simulation 
will accurately predict the physical variables of interest very well, for any of these types of 
input. As a final check, a step bump case was run and Figures 5-11 (a and b) show the 
results. It can be seen from Figure 5-1 la that the strut stroke is predicted within about 
8% of the test data. However, the model does not correctly predict when the strut locks 
up, as seen by the overshoot at around 1.25 seconds. This suggests that the criteria for 
the strut to lock up must be improved. Otherwise, the simulation predicts the response 
very well. 

Figure 5-1 lb shows the response of the upper mass to the step bump. The 
simulation predicted to within about 10% the values as recorded in the test data. The 
oscillation that can be seen in the predicted data when the strut is locked is due to tire 
bouncing. Since the test data does not show this behavior, the tire damping needs to be 
increased. 


55 




Position (in) 


Stroke Remaining vs. Time 


■ -Measured Stroke 
. Predicted Stroke 
:lnput Displacement 



Time (sec) 


Figure 5-lla: Time history of strut position as gear encounters a step bump. 


Wing/Gear Position vs. Time 




/ ^ . 


Measured Wing/Gear Positioi 
Predicted VVlng/Gear Positidr 
Input Displacement 


0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

Time (sec) 

Figure 5-llb: Time history of Wing/Gear Position as gear encounters a step bump. 
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This section has shown that the parameters found by tuning the model to a 1.0 
inch amplitude sine sweep from 0.75 to 3.75 Hz. in 40 seconds were also good for other 
cases, validating the model for other, more general types of cases that may be run. 

5.5 Summary 

In this chapter, the static parameters of the landing gear model such as system 
mass, tire load-deflection curve and pressure-stroke curve, were found from various tests 
performed at NASA Langley's Aircraft Landing Dynamics Facility. After the model had 
been updated with static data, a sine sweep test was run and the parameters of the model 
were, once again, updated to reflect the new knowledge of the system. Other tests were 
run as a check of this set of parameters and it was found that the simulation predicted the 
system response within the 10% range. Further changes of the model were suggested as 
a way to further increase the accuracy of the model. In summary, then, a fully nonlinear 
analytical model of a telescoping landing gear is presented here which has been tuned 
with test data and validated with data and which has run times of about 3 minutes per 
dynamic simulation second. For semi-static runs, times are much shorter, about 30 
seconds real time per simulated second. This model may now be adjusted for further 
accuracy or it may be used as is as a tool for evaluating the effect of applying various 
active control schemes to the landing gear. 
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Chapter 6: Concluding Remarks 


6.1 Conclusions 

The research presented in this document brings together in one place a 
comprehensive development of the equations of motion, a discussion of the problems 
associated with integrating the equations of motion, and describes a procedure to tune the 
analytical model with test data. It was found that an implicit predictor-corrector routine 
was very efficient, but a Runga-Kutta routine was needed to integrate across 
discontinuities. The result of this research is a simulation tool of the A-6 Intruder main 
landing gear which has been validated with both static and dynamic data. This model is a 
contribution to the effort to study and correct vibrations that are transmitted through the 
landing gear into the fuselage of a plane. This tool is to be used to simulate and evaluate 
the effect of various active control schemes to reduce force or vibration transmission to 
the fuselage. Simulations based on the development presented in this document currently 
exist in the form of a FORTRAN program, a SIMULINK program and in a DADS 
format. In conclusion, then, this simulation is a powerful tool that is the result of both 
analytical and experimental efforts. 

6.2 Future Research 

This research has provided a tool which can be put to many uses. A couple of 
areas of future research are suggested. The first is to complete the update of the model. 
The model stmt damping could be further tuned. A suggestion for this fine tuning is to 
run sensitivity studies on the variables that directly affect the model damping. An 
additional result may be that the friction associated with the guide roller bearings may 
need to be simulated. This tuning is optional, though, as acceptable measures of 
important parameters are currently being predicted. In addition, the stick-slip friction 
model needs to be tuned. Not only do the criteria of when to stick and when to slip need 
to be further defined, but the model that takes the dynamics from motion to the stuck 
phase and back again needs to be smoothed. With this change, the Runga-Kutta may not 
be needed at all, leaving only the Adams-Moulton, which has demonstrated excellent 
efficiency in solving the continuous portions of the total solution history. 
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The second area of future research involves actively controlled landing gear. As 
stated in the introduction, the ultimate goal of this research is to contribute to the process 
of alleviating vibrations in the HSCT. As a step toward that, this simulation is to be used 
to evaluate various control schemes. Three approaches to control aircraft vibrations using 
the landing gear are proposed. One is an active orifice concept which would allow the 
landing gear damping characteristics to be greatly altered in response to unacceptable 
ground disturbances. The second is a force actuator placed parallel to the main landing 
gear cylinder that acts to actively attenuate relative motion of the fuselage to the ground. 

A third concept uses electro-rheological fluids to actively control the damping 
characteristics of the landing gear. This approach may call for the replacement of the 
current fluid in the cylinder or perhaps be implemented as a parallel damper. The 
dynamics of these various concepts are to be modeled and added to the validated 
simulation. Control laws will then be evaluated on the basis of reducing vibrations in a 
simulated cockpit. These concepts are being developed by NASA engineers and are to be 
tested at a later date. 
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Appendix A 


A.l Summary of Program 

As part of this reasearch, a computer program, "gearfin.f ' has been written to 
numerically solve the equations of motion of landing gear in time. The notation used 
throughout the appendix is consistent with that used in the FORTRAN program, not 
that which is presented in the paper. Four states are used for the solution. These states 
are Y(l) and Y(2), the position and velocity, respectively of the wing/gear interface, and 
Y(3) and Y(4), the position and velocity, respectively, of the wheel axle. The program 
uses the Adams/Moulton predictor/corrector numerical method, and a variable step 
Runga-Kutta to get past discontinuities. 

The program's setup is simple. There are four data files associated with it. The 
first, "pin.dat", lets the user define the shape of the metering pin. The first entry is the 
number of slope changes of the pin (n). The next n entries are the n diameters. Then the 
n lengths associated with the diameters. Lastly in this data file is the number defining 
the maximum stroke. 

The next data file is called "piston.dat". This file is the entry point of the many 
(12) piston/cylinder associated parameters. The inputs of this file are: Xsi, the initial 
length at which the gear is charged, Pi, the initial charge pressure, y, the polytropic gas 
constant, Du, the diameter of the upper chamber, DL, the diameter of the lower 
chamber, D1RC and DIRE, the diameters of one snubber hole under compression and 
extension conditions respectively, Dop, the diameter of the hole in the orifice plate, and 
Dpis, the diameter of the piston shaft, Mu, the mass of the upper system, Ml, and 
FSMAX, the value of the maximum sticking friction of the piston seals. 

The third file is called "ic.dat 1 '. This is a file of the initial conditions. Inputs are 
the initial conditions of the four states in the order described above. 

The last data file is a runway profile. Since various runways could be used, the 
name of the file is left up to the user. The user needs to locate the "Read in Runway" 
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section of this program and change the name. The first line of the input file should 
contain the length (I) of the following data. The next element of the input file is the 
vector: [TIME(I) ELEV(I) ELEVD(I)]. This file contains a time history of the height of 
the runway in meters and piecewise continuous time rate of change of the height as the 
wheel rolls along it. Therefore, any type of velocity or acceleration case may be 
investigated. 

This program has six subroutines associated with it. The first is the integration 
routine "ddriv2.f '. The next is a routing called "F". Its function is to define the 
derivative of the states at any time t. The third routine is called "METPIN". This 
subroutine returns the current value of Dpin based on the current stroke and the input 
data of pin.dat. The fourth is called "COEFF". In this subroutine, the C and K 
coefficients of the landing gear equations are calculated using the input data from 
"piston.dat", Dpin, and the stroke rate. The fifth is "RKF4.P, a fourth-order Runge- 
Kutta variable step integration routine. When the program hits a discontinuity (as when 
friction suddenly changes sign) it will automatically switch from the predictor corrector 
(ddriv2.f) to Runge-Kutta (rkf4.f). The sixth routine is really a copy of "F", and is 
called "FOUT". This routine is independent of the integration and can be called to record 
the states by the user. The current data files are for the A-6 Intruder. 

The current form of the program has many outputs. These outputs are selected 
with the idea of validation of the model in mind. All output files are written in ASCII, 
and have the form of an nX4 matrix, where n is the number of times along the interval 
that the program writes the solutions. The first is "y.out". It contains the four states, 
Y(l), Y(2), Y(3), Y(4). The second file, "tfaa.out", records the time, the value of friction 
at each time, the accelerations, YDOT(2), and YDOT(4), of the two degrees of freedom. 
The third file is used to check some numerical fitness. It is called "check-out". This file 
writes the peak sticktion friction, the two forces that are subtracted to compare against 
the stiction friction, and the relative velocity of the upper and lower masses. The 
fourth file, "hydr.out", records the spring and damping coefficients of the functional 
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equations of motion. They are Pu, PI, Ps, and Qo, the upper pressure, lower presser, 
snubber pressure, and flowrate through the main orifice respectively. The fifth output 
file, "tire.out", is used to record input information and tire coefficients. These values are 
U, UDOT, K,, and C,. These output files are used in conjunction with a MATLAB m- 
file "gdplta.m" and "resp.m" to calculate and plot many variables that are of interest. 

The benefit of writing the program in a component form is that it makes the 
program much more flexible. Conceivably, many different gear configurations could be 
simulated, as long as they have the same basic setup. Also, this form of expressing the 
parameters and equations is conducive to optimization studies of various parameters. It 
is thought that studies will be performed to optimize the metering pin, perhaps the 
snubber orifices and possibly other parameters. 

This program was originally written to simulate touchdown, or drop test 
conditions. However, it was expanded to include a rollout or taxi case in which a 
runway profile is defined, with bumps at various locations. 


A.2 Program Listing 

PROGRAM GEARFIN 

C This program was developed as partial satisfaction of a 
C Masters of Science from George Washington University, Joint 
C Institute for the Advancement of Flight Sciences. Work was 
C done by James Daniels. 5/15/96. 

C 

C A computer program, "gearfin.f ' has been written to numerically solve 
C the equations of motion of landing gear in time. Four states are used 

C for the solution. These states are Y(l) and Y(2), the position and 

C velocity, respectively of the wing/gear interface, and Y(3) and Y(4), 

C the position and velocity, respectively, of the wheel axle. The 
C program uses the Adams/Moulton predictor/corrector numerical method, 
C and a variable step Runga-Kutta to get past problem discontinuities. 

C 

C The program's setup is simple. There are four data files associated 
C with it. The first, "pin.dat", lets the user define the shape of the 
C metering pin. The first entry is the number of slope changes of the 
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C pin (n). The next n entries are the n diameters. Then the n lengths 

C associated with the diameters. Lastly in this data file is the number 

C defining the maximum stroke. 

C 

C The next data file is called "piston.dat". This file is the entry 

C point of the many piston/cylinder associated parameters. The inputs of 

C this file are: Xsi, the initial length at which the gear is charged, 

C Pi, the initial charge pressure, gamma,the polytropic gas constant, 

C Du, the diameter of the upper chamber, DL, the diameter of the lower 
C chamber, D IRC and DIRE, the diameters of one snubber hole under 

C compression and extension conditions respectively, Dop, the diameter 
C of the hole in the orifice plate, and Dpis, the diameter of the piston 
C shaft, Mu, the mass of the upper system, Ml, the lower mass, DF, the 

C percent of the maximum friction that is active dynamically, and FSMAX, 

C the value of the maximum sticking friction of the piston seals. 

C 

C The third file is called "ic.dat". This is a file of the initial 

C conditions. Its inputs are the initial conditions of the four states 

C in the order described above. 

C 

C The last data file is a runway profile. Since various runways could be 
C used. The name of the file is left up to the user. The user needs to 

C locate the "Read in Runway" section of this program and change the 

C name. Its inputs are: TIME(I), ELEV(I) and ELEVD(I). This file 
C contains a time history of the height of the runway in meters as the 
C wheel rolls along it. Therefore, any type of velocity or acceleration 
C case may be investigated. 

C 

C This program has six subroutines associated with it. The first is 

C the integration routine "ddriv2.f '. The next is a routing called "F". 

C Its function is to define the derivative of the states at any time t. 

C The third routine is called "METPIN". Its job is to return the current 
C value of Dpin based on the current stroke and the input data of 
C pin.dat. The fourth is called "COEFF". It calculates the C and K 
C coefficients of the landing gear equations using the input data from 
C "piston.dat", Dpin, and the stroke rate. The fifth is "RKF4.f', a 
C fourth-order Runge-Kutta variable step integration routine. 

C When the program hits a discontinuous spot (as when friction suddenly 
C changes sign) it will automatically switch from the predictor 
C corrector (ddriv2.f) to Runge-Kutta (rkf4.f). The sixth routine is 
C really a copy of "F", and is called "FOUT". This routine is independent 
C of the integration and can be called to record the states by the user. 

C 
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c 


The current data files are for the A- 6 Intruder. 


EXTERNAL F,COEFF,METPIN 


C Variable Declaration 

C — Main — 

DOUBLE PRECISION Y, YDOT, WORK,DT 
DOUBLE PRECISION T, TOUT, EWT, EPS 
DOUBLE PRECISION G,STP,TMIN 
INTEGER NEQ,L,I,J,JJ 

INTEGER MSTATE, MINT, LENW,IWORK, LOUT 
INTEGER LENIW,N,NROOT,COUNT,NUMBR 

C — Subroutine F and FOUT — 

DOUBLE PRECISION CT,KT,LIFT,TEMP 
DOUBLE PRECISION MU,GRAV,ML,U,UDOT,FRICT 
DOUBLE PRECISION XS,FT,FC,MA 
DOUBLE PRECISION DELTA, FR,MTIRE,FTC,FTK 
DOUBLE PRECISION F1,F2,FSTICK,DEL 

C —Define Type for Input Calculation— 

DOUBLE PRECISION PIE, ELVU,ELVL,ELEV, TIME 
DOUBLE PRECISION DISTL,DISTU,TAXISPD,DIS 
DOUBLE PRECISION AA,BB,HGT,LE,AO,DD,X,ELEVD 
INTEGER LN 

C — Subroutine METPIN — 

DOUBLE PRECISION PAR,DPIN,D,LNG 
INTEGER NUM 

C — Subroutine COEFF — 

DOUBLE PRECISION PARl,C,K,CON,ALS 
DOUBLE PRECISION DOR,AO,RHO,CD,AL,FLW 
DOUBLE PRECISION AOP,APIN,MEW,VEL,RD,BETA,C 1 
DOUBLE PRECISION ARC,ARE,CDE,CDC,E1,E2,E3,E4,AR 
DOUBLE PRECISION BETAE,BETAC,AS 

C — Subroutine RKF4 — 

DOUBLE PRECISION TOL,PD,HMIN,HMAX,H,WK 
INTEGER MTH,IERR 
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C — Define Array Size — 

PARAMETER (NEQ=4) 

DIMENSION Y (NEQ), YDOT (NEQ), WK(7*NEQ) 

DIMENSION W ORK(3 00),I W ORK(30) 

DIMENSION P AR1 ( 1 2),C(2),K(2) 

DIMENSION PAR(20),D( 1 0),LNG( 1 0),FLW(4) 

DIMENSION TIME(5000),ELEV(5000),ELEVD(5000) 

COMMON /P ARAM/P AR 1 , PAR, STP, TIME, ELEV,ELEVD,NUM,LN 

C Read in Metering Pin Info 

OPEN(UNIT=4, STATUS- OLD', FILE-pin.dat') 

READ(4,*) NUM 
DO I=1,2*NUM 
READ(4,*) PAR(I) 

END DO 

READ(4,*) PAR(2*NUM+1) 

CLOSE(UNIT=4, STATUS- KEEP') 

C P AR(2 *NUM+ 1 )=[ D(N), L(N), XSMAX ] 

C Read in Piston Info for Coefficients 

OPEN(UNIT=5 , ST ATUS- OLD', FILE- piston.dat') 

DO 1=1,12 
READ(5,*) PAR 1(1) 

END DO 

CLOSE(UNIT=5 ,ST ATU S- KEEP') 

C PAR1(13)=[ XSI,PI, YI,DU,DL,D 1 RC,D 1 RE,DOP,DPIS MU ML FSMAX ] 

C Read in Runway Profile 

C The data file is an ascii, three column file that 

C contains time in the first column and runway height 
C in meters in the second column and time derivative of 
C runway height in third column. 

OPEN(UNIT=6, STATUS-UNKNOWN', FILE- r2 1 1 .dat’) 

READ(6,*) LN 
DO 1=1, LN 

READ(6,*) TIME(I), ELEV(I), ELEVD(I) 

END DO 

CLOSE(UNIT=6, STATUS-KEEP') 


C- 


-Define Height of Stop Block in cylinder — 
STP=9.5/39.37 



O U 


C Open Output Files 

C These variables are strictly up to the user. These 
C are the current variables being recorded. A MATLAB 
C routine "gdplta.m" exists to plot these particular 
C variables in a coherent way. It may be altered by 
C any user. 

C WRITE(1 1,790) (Y(J),J=1,4) 

OPEN(UNIT=l 1 , S T ATU S-UNKN O WN',FILE='y .out’) 

C WRITE(1 2,790) T,FRICT,YDOT(2),YDOT(4) 

OPEN(UNIT= 1 2, STATUS- UNKN O WN',FILE='tfaa.out') 

C WRITE (13,790) FR, FI, F2,VEL 

OPEN(UNIT= 1 3, STATUS='UNKNOWN', FILE- check.out') 
C WRITE(14,790) (FLW(J),J=1,4) 

OPEN(UNIT= 1 4, ST ATU S-UNKN OWN',FILE=’hydr.out') 
C WRITE(15,790) U,UDOT,KT,CT 

OPEN(UNIT= 1 5 , STATU S- UNKNO WN',FILE- tire.out') 

C 

C — Landing or AM-2 Mat bump case or runway profile — 

LOUT= 160000 
DT=0. 00025 

C — The amount of time one gets is = DT*LOUT 


Re-initialize the initial conditions. 

DO 1=1, NEQ 
Y(I)=0.0 
YDOT(I)=0.0 
END DO 

C Read initial Conditions 

C Two files are available, "landic.dat" or "static.dat". The first 
C startswith positional vector at zero and velocity vector at a given 
C sink rate. The second starts with the position vector in its 
C equilibrium position and the velocity vector set to zero. 

OPEN(UNrr=10,STATUS-OLD',FILE- r211ic.dat') 
READ(10,*) (Y(I),I=1,4) 

CLOSE(UNIT= 1 0, ST ATUS-KEEP’) 

T=0.0 



nnn 


C Call into current memory current values 
C of all the states and their accelerations. 

CALL F(N,T,Y,YDOT) 

100 TOUT=0.0 

C Define DDRIV2 Parameters 

C See the leading introduction of DDRIV2.f for further 

C explanation of these and other parameters. 
MSTATE=1 
NROOT=0 
EPS=lE-6 
EWT=1E-15 
MINT=3 
LENW=300 
LENIW=30 
COUNT=0 


Loop entire process to get LOUT data points. 


DO L=l,LOUT 

C Increment time step. 

TOUT=DT*L 

C Call main integration routine. 

5 CALL DDRIV2(N,T,Y,F, TOUT, MSTATE, NROOT,EPS,EWT, MINT, WORK, 
* LENW,IWORK,LENIW,G) 


C Provide a visual check of the progress of the integration. 

IF ((REAL(L)/2000.0) .EQ. INT(REAL(L)/2000.0)) THEN 
WRITE (*,889) L 
END IF 

C — Check for errors from DDRIV2. 

C Error (3) or (-3) indicates too much work to ge to next 
C time step. Time to switch to Runga-Kutta. 

IF (MSTATE .EQ. 3 .OR. MSTATE .EQ. -3) THEN 

C Initiate backup integration method since DDRTV2 cannot 
C complete the next integration step. 


70 



C RKF4 Parameters 

T0L=lE-5 

MTH=1 

HMIN=1E-15 

HMAX=lE-2 

H=HMAX 

CALL F(N,T,Y,YDOT) 

2 CALL 

RKF4(NEQ,T, TOUT, Y,TOL,F, PD, MTH,HMIN,HMAX,H,WK, IERR) 
C Check for error from R-K. (-1) indicates that the tolerances 
C are not set up correctly for this problem. 

IF (IERR .EQ. -1) THEN 
WRITE (*,889) IERR 
TOL=lE-5 
MTH=1 
HMIN=1E-15 
HMAX=lE-2 
H=HMAX 
IERR=0 

IF (COUNT .EQ. 500) GOTO 500 
COUNTCOUNT+1 
GOTO 2 
ENDIF 

IF (IERR .EQ. 0) THEN 
MSTATE=1 
GOTO 1 
ENDIF 
ENDIF 

1 COUNT=0 


C — Define some filtering process to record data. Otherwise, data 
C — files are VERY large. 

IF ((RE AL(T /DT)/40 . 0) .EQ. INT(REAL(T/DT)/40.0)) THEN 
CALL FOUT(N,T,Y,YDOT) 

END IF 

C END THE L LOOP: i.e. Each integration step 
END DO 

889 FORMAT(I6) 

500 CLOSE(UNIT=l 1, STATUS- KEEP') 

CLOSE(UNIT=12, STATUS- KEEP') 
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CLOSE(UNIT= 1 3, STATUS- KEEP') 
CLOSE(UNIT = 1 4,ST ATU S- KEEP') 
CLOSE(UNIT = 1 5,STATUS- KEEP') 


C— End of Main program. 
END 


C 

C=== 

c 

C Subroutines 

C 

C— — — — — 

c 


SUBROUTINE F(N,T,Y,YDOT) 


C Subroutine F defines the functional form of Ydot(i)= ... 

C Anything changed in this subroutine needs also to be changed 
C in the FOUT routine and vise versa. 

EXTERNAL COEFF,METPIN 

DOUBLE PRECISION Y,YDOT,T,CT,KT,LIFT,MU,GRAV 
DOUBLE PRECISION ML,U,UDOT,FRICT,DPIN,TEMP 
DOUBLE PRECISION XS,VEL,FT,FC,MA,DELTA,FR 
DOUBLE PRECISION MTIRE,DEL,STP,F 1 ,F2,FTC,FTK 
DOUBLE PRECISION FSTICK,PAR,PAR1,C,K,FLW 
C —Define Type for Input Calculation— 

DOUBLE PRECISION PIE,ELVU,ELVL,ELEV,TIME 
DOUBLE PRECISION DISTL,DISTU,TAXISPD,DIS 
DOUBLE PRECISION AA,BB,HGT,LE,AO,DD,X,ELEVD 

INTEGER N,NUM,LN 

PARAMETER (NEQ=4) 

DIMENSION Y(NEQ),YDOT(NEQ),TIME(5000),ELEV(5000),ELEVD(5000) 
DIMENSION C(2),K(2),PAR(20),PAR 1 ( 1 2),FLW(4) 
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COMMON /P ARAM/P AR1 , PAR, STP,TIME,ELEV,ELEVD,NUM,LN 


MU=PAR1(10) 

ML=PAR1(1 1) 

XS=P AR(2 *NUM+ 1 )+Y (3 )- Y ( 1 ) 

CALL METPIN(Y,NUM,PAR,DPIN) 

CALL COEFF(Y,PARl,DPIN,XS,C,K,FLW) 


C — Calculate the input U into the tire from runway — 

C This section defines the ground excitation for a runway profile 
C input case. It linearly interpolates between points of the data file. 

DO I=1,LN-1 

IF(T .GE. TIME(I) .AND. T .LT. TIME(I+1)) THEN 
ELVL=ELEV(I) 

EL VU=ELE V (1+ 1 ) 

ELLD=ELEVD(I) 

ELUD=ELE VD(I+ 1 ) 

U=(ELVU-ELVL)*(T-TIME(I))/(TIME(I+ 1 )-TIME(I))+ELVL 
C UDOT is dU/dT, leaving only: 

UDOT=(ELUD-ELLD)*(T-TIME(I))/(TIME(I+l)-TIME(I))+ELLD 
ENDIF 
END DO 
C 

C — Define Input of AM2 Repair Mat — 

C TAXISPD=20.0*.5 144444 

C AA= 15.0*(.3048) 

C BB= AA+4.0*(.3048) 

C CC= BB+70.0*(.3048) 

C DD= CC+4.0*(.3048) 

C HGT= 1.5/39.37 

C - X=TAXISPD*T - 

C X=TAXISPD*T 

C IF (X .LT. AA .OR. X .GT. DD) THEN 

C U=0.0 

C UDOT=0.0 

C ELSEIF (X .GE. AA .AND. X .LE. BB) THEN 

C U=HGT*(X-AA)/(BB-AA) 

C UDOT=HGT*TAXISPD/(BB-AA) 
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C ELSEIF (X .GT. BB .AND. X .LT. CC) THEN 
C U=HGT 

C UDOT=0.0 

C ELSEIF (X .GE. CC .AND. X .LE. DD) THEN 
C U=-HGT*(X-CC)/(DD-CC)+HGT 

C UDOT=-HGT*TAXISPD/(DD-CC) 

C END IF 

C — Toggle for landing/runway input case. We want no input 
C — for landing case, but DO want bumps etc. for other cases. 

C U=0.0 

C UDOT=0.0 

C Tire Model as updated from Experimental Data — 

KT=(-252.0*(39.37*(Y(3)+U))**2.0+1397.0*(Y(3)+U)*39.37+4267.0) 
* *(39.37*4.4482) 

C — Tire Damping model, as observed from test data — 

MTIRE=ML 

CT=5000.0 

C — Define the Tire Force (FT) — 


FTK=1.0*KT*(Y(3)+U) + 130.0*4.4482 
FTC=CT*(1 ,0)*(Y(4)+UDOT) 

IF ((Y(3)+U) .LT. 0.0) THEN 
FTK=0.0 
FTC=0.0 
END IF 
FT=FTC+FTK 
C 


GRAV=9.81 


— Lift Model — 

LIFT=9.81*MU 

LIFT=0.0 

VEL=Y(2)-Y(4) 


C Defining relative forces before friction — 

F 1 =MU*GRAV-LIFT+C( 1 )* VEL* *2.0+K( 1 )*( 1 .0/XS)**(PARl (3)) 
F2=ML*GRAV+C(2)*VEL**2.0+K(2)*(1.0/XS)**(PAR1(3))-FT 


C 

C Add the KARNOPP friction Model to the accelerations. 
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C —DEL is how close relative velocity needs to be to zero to stick. 

DEL=.0009 
DELTA=ABS(F 1 -F2) 

C Calculate the bearing friction of the piston in the 
C cylinder. FC is frictional coefficient, MA is moment arm. 

C FC = .05 

C MA= 10.5/39.37 

C TEMP=(F C *FT* MA/(AB S( Y ( 1 )- Y(3 )+STP))+P AR 1(12))*. 75 

C — Future Friction Model. Needs to be ironed out. 

C TEMP=(4000.0*EXP(-XS/( 1 .0* .0254))+ 1 000.0* 

C * EXP(-ABS(VEL/0.05)))*4.44822 
C * *EXP(-ABS(VEL/0.05)) 

C FRICT=-TANH(VEL/.008)*TEMP 

C Also part of future friction model. 

C FR=1.0*(4000.0*EXP(-XS/(6.2*.0254))+1000.0*.4)*4.4482+PARl(12) 

C Friction toggle for finding initial conditions — 

TEMP=400.0 

FRICT=(-0.0-TANH(VEL/.008))*TEMP*4.4482 
FR=2689. 0*4.44822 
C FR=0.0 

C FRICT=0.0 

C 


IF (DELTA .LT. FR .AND. ABS(VEL) .LT. DEL) THEN 
C Case 1, Piston Sticks in Cylinder. 

F1=MU*GRAV-LIFT+K(1)*(1.0/XS)**(PAR1(3)) 
F2=ML * GRA V +K(2) * (1.0/XS)**(PAR1 (3 ))-FT 
FSTICK=(ML*F 1 - MU*F2)/(MU+ML) 

YDOT(l)=Y(2) 

YDOT (2)=F 1 /MU - FSTICK/MU 
YDOT(3)=Y(4) 

YDOT(4)=F2/ML + FSTICK/ML 
ELSE 

C Case 2, Relative Motion between Piston 
C and Cylinder, with friction present. 

YDOT(l)=Y(2) 

YDOT (2)=F 1 /MU + FRICT/MU 
YDOT(3)=Y(4) 

YDOT(4)=F2/ML - FRICT/ML 
END IF 
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C Case 3, The gear and tire leave the ground. 

IF ((Y(3)+U) .LT. 0.0 .AND. XS .GT. PAR(2*NUM+1)) THEN 
YDOT(l)=Y(2) 

YD0T(2)=F1/MU - FSTICK/MU 
YDOT(3)=Y(2) 

YDOT (4)=F2/ML + FSTICK/ML 
END IF 

RETURN 

STOP 

END 


C SUBROUTINE METPIN 

SUBROUTINE METPIN(Y,NUM,PAR,DPIN) 

C This subroutine linearly interpolates in the D(n), diameter 
C array, and L(n), the length array to determine the diameter 
C of the metering pin, DPIN, at any stroke, XS. 

DOUBLE PRECISION Y,PAR,DPIN,D,LNG,XS 
INTEGER NUM.I 
PARAMETER (NEQ=4) 

DIMENSION Y (NEQ),PAR(20),D( 1 0),LNG( 1 0) 

C PAR(2*N+1) = [ D(N),LNG(N),XSMAX ] 

DO 1=1, NUM 
D(I)=PAR(I) 

LN G(I)=P AR(NUM+I) 

END DO 

DPIN=D(NUM) 

XS=P AR(2 *NUM+ 1 )- Y( 1 )+ Y(3) 

DO I=1,NUM-1 

IF (XS .LT. LNG(I) .AND. XS .GT. LNG(I+1» THEN 
DPIN = D(I) + (D(I+1)-D(I))*(LNG(I)-XS)/(LNG(I)-LNG(I+1)) 
ENDIF 
END DO 
RETURN 

END 

C SUBROUTINE COEFF 
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SUBROUTINE COEFF(Y,PARl ,DPIN,XS,C,K,FLW) 


C For a thourough understanding of this subroutine, look 

C in the thesis Chapter 2, Equations (1) and (2) as they 
C are presented at the end of the chapter, this subroutine 
C calculates the coefficients as they are defined in the paper. 

C These coefficients are only functions of geometry and flow 
C direction. 

DOUBLE PRECISION Y,PARl,C,K,CON,ALS 

DOUBLE PRECISION DOR,AO,RHO,CD,AL,FLW 

DOUBLE PRECISION DPIN,AOP,APIN,XS 

DOUBLE PRECISION MEW,VEL,RD,BETA,C 1 

DOUBLE PRECISION ARC,ARE,CDE,CDC,E1,E2,E3,E4,AR 

DOUBLE PRECISION BETAE,BETAC,AS 

INTEGER N 

PARAMETER (NEQ=4) 

DIMENSION Y(NEQ),PAR1( 1 2),C(2),K(2),FLW(4) 

C PAR1(9)=[ XSI,PI, YI,DU,DL,D 1 RC,D 1 RE,DOP,DPIS,MU,ML FSMAX ] 

C — Calculate Various areas to be used — 

CON=7853981 
ARC=CON* P AR 1 (6) * *2 .0 
ARE=CON*PAR1(7)**2.0 
AL=CON*PAR1(5)**2.0 
AS=AL-CON*PAR1(9)**2.0 
AOP=CON*PAR1(8)**2.0 
APIN = CON*DPIN**2.0 
AO=AOP-APIN 


C Define a constant (Cl) to account for the annular 
C nature of DOR. It is effectively SMALLER to 
C the fluid. I'm not using this at the moment, but 
C should the tests bare out the fact that something 
C is not working quite right, this may be a parameter 
C to tweak. 

Cl=1.0 

DOR=C 1 *SQRT(AO/.785398 1 ) 

C RHO is fluid density. MEW is fluid viscosity. 

RHO=912.0 


77 



MEW=35.0*.001 
VEL=Y (2)- Y (4) 

C RD is reynolds number. An RD model of the discharge 
C coefficient may want to be used in the future. 

RD=RHO*ABS(VEL)*PARl(5)/MEW 
C The various BETA'S are ratios of (fluid coming from Dl)/ 

C (fluid going through D2) — > (Dl/D2)=Beta 
BETA=DOR/PARl(5) 

BETAC=P AR 1 (6)/P AR 1 (5) 
BETAE=PAR1(7)/(SQRT(AS/(12.0*.7853981))) 

C Therefore, the discharge coefficients are now only functions 
C of geometry, not Reynold's number. 

CD=1.0*BETA**2.0 - .4813*BETA + .8448 
CDC=.95*(.8*BETAC**2.0 - .4813*BETAC + .8448) 
CDE=1.0*(.8*BETAE**2.0 - ,4813*BETAE + .8448) 

C COMPRESSION 

IF (VEL .GT. 0.0 .OR. VEL .EQ. 0.0) THEN 
C AR is the sum of the twelve areas that comprise the snubber 
C orifices. 

AR=12.0*ARC 

E1=AO*.95*CD*SQRT((2.0)/(RHO*(1.0-BETA**4.0))) 
E2=AR*CDC*SQRT(2.0/(RHO*( 1 ,0-BETAC**4.0))) 

ALS=AL-AS 

C(1)=((ALS/El)**2.0-(AS/E2)**2.0)*AS-(ALS/E1)**2.0*(AL-A0) 
C(2)=((AS/E2)**2.0-(ALS/E1)**2.0)*(AS-AR)+(ALS/E1)**2.0*(AL-AR) 
K( 1 )=(AS-AL)*PAR 1 (2)*PAR 1 ( 1 )**PAR1 (3) 

K(2)=ALS*PAR 1 (2) *P AR 1(1)* *P AR 1 (3) 

C~As a bonus, the pressures and flow rates are also being calculated here. 
PU=PAR1(2)*(PAR1(1)/XS)**PAR1(3) 

PL=PU+(((AL-AS)/E 1 )*' VEL)* *2.0 
PS=PL-(AS*VEL/E2)**2.0 
QO—E1 *SQRT(PL-PU) 

C EXTENSION 
ELSE 

C AR is the sum of the twelve areas that comprise the snubber 
C orifices. 

AR=12.0*ARE 

E3=AO* 1 . 1 *CD*SQRT(2.0/(RHO*( 1 .0-BETA* *4.0))) 
E4=AR*CDE*SQRT(2.0/(RHO*(1.0-BETAE**4.0))) 

ALS=AL-AS 

C(1)=(ALS/E3)**2.0*(AL-A0)+((AS/E4)**2.0-(ALS/E3)**2.0)*AS 

C(2)=((ALS/E3)**2.0-(AS/E4)**2.0)*(AS-AR)-(ALS/E3)**2.0*(AL-AR) 
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K( 1 )=(AS-AL)*PAR 1 (2)*PAR 1 ( 1 )* *PAR1 (3) 

K(2)=ALS * PAR 1 (2) *P AR 1 ( 1 ) * * P AR 1 (3) 

C— As a bonus, the pressures and flow rates are also being calculated here. 
PU=P AR 1 (2) *(PAR 1(1 )/XS) * * PAR 1(3) 

PL=PU-((( AL- AS)/E3) * VEL) **2.0 
PS=PL+(AS*VEL/E4)**2.0 
QO=E3 * SQRT (PU-PL) 

END IF 

C— Put the pressures and flow rates into an array to pass 
C~outward for recording. 

FLW(1)=PU 

FLW(2)=PL 

FLW(3)=PS 

FLW(4)=QO 

RETURN 

END 

C 


SUBROUTINE FOUT(N,T,Y,YDOT) 


C Subroutine FOUT is the same as F and comments match until the output. 
C Anything changed in this subroutine needs also to be changed 
C in the F routine and vise versa. 

EXTERNAL COEFF.METPIN 

DOUBLE PRECISION Y,YDOT,T,CT,KT,LIFT,MU,GRAV 
DOUBLE PRECISION ML,U,UDOT,FRICT,DPIN,TEMP 
DOUBLE PRECISION X S, VEL,FT,FC,MA,DELTA,FR 
DOUBLE PRECISION MTIRE,DEL,STP,F 1 ,F2,FTC,FTK 
DOUBLE PRECISION FSTICK,PAR,PAR1,C,K,FLW 
C —Define Type for Input Calculation— 

DOUBLE PRECISION PIE,ELVU,ELVL,ELEV,TIME 
DOUBLE PRECISION DISTL,DISTU,TAXISPD,DIS 
DOUBLE PRECISION AA,BB,HGT,LE,AO,DD,X,ELEVD 

INTEGER N,NUM,LN 

PARAMETER (NEQ=4) 
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DIMENSION Y(NEQ),YDOT(NEQ),TIME(5000),ELEV(5000),ELEVD(5000) 
DIMENSION C(2),K(2),PAR(20),PAR 1 ( 1 2),FLW(4) 

COMMON /P ARAM/PAR 1,PAR,STP, TIME, ELEV,ELEVD,NUM,LN 

MU=PAR1(10) 

ML=PAR1(1 1) 

XS=P AR(2*NUM+ 1 )+' Y(3 )-' Y( 1 ) 

CALL METPIN(Y,NUM,PAR,DPIN) 

CALL COEFF(Y,PARl ,DPIN,XS,C,K,FLW) 


C — Calculate the input U into the tire from runway — 

C This section defines the ground excitation for a runway profile 
C input case. It linearly interpolates between points of the data file. 

DO I=1,LN-1 

IF(T .GE. TIME(I) .AND. T .LT. TIME(I+1)) THEN 
ELVL=ELEV(I) 

ELVU=ELEV(I+1) 

ELLD=ELEVD(I) 

ELUD=ELE VD(I+ 1 ) 

U=(ELVU-ELVL)*(T-TIME(I))/(TIME(I+ 1 )-TIME(I))+ELVL 
UDOT=(ELUD-ELLD)*(T-TIME(I))/(TIME(I+l)-TIME(I))+ELLD 
ENDIF 
END DO 
C 

C — Define Input of AM2 Repair Mat — 

C TAXISPD=20.0*.5 144444 

C AA= 15.0*(.3048) 

C BB= AA+4.0%3048) 

C CC= BB+70.0%3048) 

C DD= CC+4.0*(.3048) 

C HGT= 1.5/39.37 

C - X=TAXISPD*T - 

C X=TAXISPD*T 

C IF (X .LT. AA .OR. X .GT. DD) THEN 

C U=0.0 

C UDOT=0.0 

C ELSEIF (X .GE. AA .AND. X .LE. BB) THEN 

C U=HGT* (X- AA)/(BB - AA) 
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C UDOT=HGT*TAXISPD/(BB-AA) 

C ELSEIF (X .GT. BB .AND. X .LT. CC) THEN 
C U=HGT 

C UDOT=0.0 

C ELSEIF (X .GE. CC .AND. X .LE. DD) THEN 
C U=-HGT*(X-CC)/(DD-CC)+HGT 

C UDOT=-HGT*TAXISPD/(DD-CC) 

C END IF 

C 

C — Toggle for landing/runway input case. We want no input 
C — for landing case, but DO want bumps etc. for other cases. 

C U=0.0 

C UDOT=0.0 

C Tire Model as updated from Experimental Data — 

KT=(-252. 0*(39.37*(Y(3)+U))**2.0+1 397. 0*(Y(3)+U)*39. 37+4267.0) 
* *(39.37*4.4482) 

C — Tire Damping model, as observed from test data — 

MTIRE=ML 

CT=5000.0 

C — Define the Tire Force (FT) — 


FTK=1.0*KT*(Y(3)+U) + 130.0*4.4482 
FTC=CT*(1 .0)*(Y(4)+UDOT) 

IF ((Y(3)+U) .LT. 0.0) THEN 
FTK=0.0 
FTC=0.0 
END IF 
FT=FTC+FTK 
C 


GRAV=9.81 

— Lift Model — 

LIFT=9.81*MU 

LIFT=0.0 


VEL=Y(2)-Y(4) 


C Defining relative forces before friction — 

F1=MU*GRAV-LIFT+C(1)*VEL**2.0+K(1)*(1.0/XS)**(PAR1(3)) 

F2=ML*GRAV+C(2)*VEL**2.0+K(2)*(1.0/XS)**(PAR1(3))-FT 


C 
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C Add the KARNOPP friction Model to the accelerations. 

C --DEL is how close relative velocity needs to be to zero to stick. 
DEL=.0009 
DELT A=AB S(F 1 -F2) 

C Calculate the bearing friction of the piston in the 
C cylinder. FC is frictional coefficient, MA is moment arm. 

C FC = .05 

C MA = 10.5/39.37 

C TEMP=(FC*FT*MA/(ABS(Y(1)-Y(3)+STP))+PAR1(12))*.75 

C — Future Friction Model. Needs to be ironed out. 

C TEMP=(4000.0*EXP(-XS/(1.0*.0254))+1 000.0* 

C * EXP(-ABS(VEL/0.05)))*4.44822 
C * *EXP(-ABS(VEL/0.05)) 

C FRICT==-TANH(VEL/.008)*TEMP 

C Also part of future friction model. 

C FR=1.0*(4000.0*EXP(-XS/(6.2*.0254))+1000.0*.4)*4.4482+PAR1(12) 

C Friction toggle for finding initial conditions — 

TEMP=400.0 

FRICT=(-0.0-TANH(VEL/.008))*TEMP*4.4482 
FR=2740. 0*4.44822 
C FR=0.0 

C FRICT=0.0 

C 


IF (DELTA .LT. FR .AND. ABS(VEL) .LT. DEL) THEN 
C Case 1 , Piston Sticks in Cylinder. 

F 1 =MU * GRA V -LIFT +K( 1 )*( 1 .0/XS)* *(PAR1 (3)) 
F2=ML*GRAV+K(2)*(1.0/XS)**(PAR1(3))-FT 
FSTICK=(ML*F1 - MU*F2)/(MU+ML) 

YDOT(l)=Y(2) 

YDOT(2)=Fl/MU - FSTICK/MU 
YDOT(3)=Y(4) 

YDOT (4)=F2/ML + FSTICK/ML 
ELSE 

C Case 2, Relative Motion between Piston 
C and Cylinder, with friction present. 

YDOT(l)=Y(2) 

YDOT (2)=F 1 /MU + FRICT/MU 
YDOT(3)=Y(4) 

YDOT (4)=F2/ML - FRICT/ML 
END IF 
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C Case 3, The gear and tire leave the ground. 

IF ((Y(3)+U) .LT. 0.0 .AND. XS .GT. PAR(2*NUM+1» THEN 
YDOT(l)=Y(2) 

YD0T(2)=F1/MU - FSTICK/MU 
YDOT(3)=Y(2) 

YDOT(4)=F2/ML + FSTICK/ML 
END IF 

WRITE(I1,790) (Y(J),J=1,4) 

WRITE( 1 2,790) T,FRICT,YDOT(2),YDOT(4) 

WRITE(13,790) FR,F1,F2,VEL 
WRITE( 14,790) (FLW(J),J=1,4) 

WRITE( 15,790) U,UDOT,KT,CT 

790 F0RMAT(E14.4,1X,E14.4,1X,E14.4,1X,E14.4) 

RETURN 

600 CLOSE(UNIT=l 1, STATUS- KEEP') 

CLOSE(UNIT=12, STATUS- KEEP') 

CLOSE(UNIT=l 3, STATUS- KEEP') 

CLOSE(UNIT=l 4, STATUS- KEEP') 

CLOSE(UNIT=l 5, STATUS- KEEP') 

STOP 

END 


A.3 Sample Input Files 

Pin.dat: 

6 

.0133604 

.021844 

.022352 

.022352 

.026162 

.026162 

0.4461 

0.353314 

0.277114 

0.112014 

0.035814 

-0.0254 

.383286 
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This file describes a metering pin of an A-6 intruder. There are six slope changes on the 
pin, and the diameters are listed from top of pin (not piston head end) to the piston end. 
Following these six diameters are the six stroke lengths associated with each diameter. 


Piston.dat: This file contains the twelve parameters as described in the program 
summary section. 

.0889 

2571744.47 

1.1 

.1524 

.1524 

3.9878 le-3 

1.587503e-3 

.0285877 

.1397 

4139.8841 

145.1 

511.5455 

Ic.dat: This file contains the initial conditions of the state vector. 

0.3164 

0.0 

0.04045 

0.0 


Test.dat: This file contains a sample runway. Only a few of the 1640 entries are shown. 


1640 

0.0000000e+00 

2.5000000e-02 

5.0000000e-02 

7.5000000e-02 

1.0000000e-01 

1.2500000e-01 

1.5000000e-01 


-2.1399094e-02 
-1. 860842 le-02 
-1. 57488 13e-02 
-1.259845 le-02 
-9.1 9330 14e-03 
-5.8510932e-03 
-2.4309578e-03 


2.7217370e-01 
2.7906726e-01 
2.8596082e-01 
3.1503620e-01 
3.405 1494e-01 
3.3422082e-01 
3.4201354e-01 


etc. 



A.4 Output Manipulation File 

% This file loads, manipulates and plots the 
% output data from the simulation. This is 
% a MATLAB .m file and is consistant with 
% MATLAB release 4.2c. 

load y.out 
load tfaa.out 
load check-out 
load hydr.out 
load tire.out 


t=tfaa(:,l); 

fr=check(:,l)*.2248089; 

fl=check(:,2); 

f2=check(:,3); 

delta=.5*abs(fl -f2)*.2248089; 

vel=check(:,4); 

xsmax=.383286; 

MU=4 139. 8841; 

ML=145.1; 

fwg=MU*(-tfaa(:,3)+9.8 1)*.2248; 

ges— tfaa(:,3); 

s=date; 

xwg=y(:,l); 

xa=y(:,3); 

vwg=y(:,2); 

va=y(:,4); 

xs=(xsmax-(xwg-xa))*39.3 7 ; 

kt=tire(:,2); 

ct=tire(:,3); 

k=hydr(:,l); 

cl=hydr(:,3); 

c2=hydr(:,4); 

frict=tfaa(:,2); 

u=tire(:,l); 

relvel=(vwg-va)*3. 28084; 

pu=hydr(:,l)*1.450377e-4; 
pl=hydr(:,2)*1.450377e-4; 
ps=hy dr( :,3)*1.450377e-4; 
qo=hy dr( : ,4) * 2 64 . 1 72052; 
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% Each of the following plots are 
% optional. They are only a sample 
% of the types of things that can be 
% considered with the information 
% available. 

subplot(2, 1 , l),plot(t,ges,'y') 
xlabel(Time (sec)') 
ylabel('Awg (g)') 
titleCWing/Gear Force vs. Time') 
gtext(s) 
grid 

subplot(2 , 1 ,2),plot(t,xs,t,u* 3 9.37,'—') 
xlabel(Time (sec)') 
ylabel('Stroke Remaining (in)') 
title('Stroke Remaining vs. Time') 
grid 

legend('-', 'Stroke','—', 'Input Displacement ') 


%figure 

%subplot(2, 1,1), plot(t,relvel,'y') 
%xlabel('Time (sec)') 

%ylabel('Relative Vel. (ft/s)’) 

%title('Relative Velocity vs. Time') 

%gtext(s) 

%grid 

%subplot(2, 1 ,2), plot(t, delta,'-', t,fr,’—') 
%xlabel(’Time (sec)') 

%ylabel('Force (lbf)') 

%title('Relative Force and Fr vs. Time') 
%grid 

%legend('-','Fwg-Fa','— ','Peak Friction Force') 

%figure 

%plot(t,pu) 

%gtext(s) 

%xlabel(Time (sec)') 

%ylabel('Pneumatic Press, (psi)') 
%titleCNitrogen Pressure vs. Time') 

%grid 
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%figure 

%plot(t,pl) 

%gtext(s) 

%xlabel('Time (sec)') 

%ylabel('Hydraulic Press, (psi)') 

%title('Fluid Pressure above Piston vs. Time') 
%grid 

%figure 

%plot(t,ps) 

%gtext(s) 

%xlabel('Time (sec)') 

%ylabel('Hydraulic Press, (psi)') 

%title(Tluid Pressure in Snubber vs. Time') 
%grid 

%figure 

%plot(t,qo) 

%gtext(s) 

%xlabel('Time (sec)') 

%ylabel('Flow Rate (gal/s)’) 

%title('Flow through Main Orifice vs. Time') 
%grid 
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